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Abstract 

The  use  of  decoys  in  combat  has  become  more  advanced  in  recent  years.  Some  of 
the  newest  military  aircraft,  such  as  the  US  Navy’s  F/A-18E/F  Superhomet,  have  the 
capability  to  deploy  a  towline  with  an  attached  decoy  when  entering  hostile  territory  as  a 
defense  mechanism  against  enemy  threats.  In  steady  state,  the  towline  extends  behind 
and  below  the  aircraft.  A  major  concern  is  the  position  of  the  towline,  as  aircraft 
maneuvers  can  cause  the  line  to  enter  the  engine  plume.  The  high  exhaust  heat  can  cause 
problems,  such  as  damaging  electrical  equipment  and  severing  the  line.  In  order  to  better 
understand  the  behavior  of  the  towline,  as  well  as  setting  up  a  method  to  analyze  the  heat 
transfer  to  the  towline,  computer  modeling  has  been  utilized  using  numerical  integration 
with  the  method  of  characteristics. 

The  method  of  characteristics  has  been  applied  to  4  hyperbolic  equations  of 
motion,  leaving  2  parabolic  equations  of  motion  to  be  calculated  at  each  timestep.  The 
energy  equation  for  heat  transfer  to  the  towline  was  also  derived,  which  provides  a  means 
to  find  local  air  density  and  towline  temperature.  From  these  a  model  was  created  to 
observe  towline  behavior  and  temperature,  which  is  shown  to  be  consistent  with  past 
research.  This  model  is  applicable  to  any  towed  body  in  any  medium  with  zero  slack 
conditions. 

The  effects  of  transient  aircraft  maneuvers  on  towline  behavior  in  a  predetermined 
temperature  field  were  analyzed  under  different  conditions  using  a  code  developed  in 
MATLAB®.  This  code  is  included  such  that  aircraft  maneuvers  in  unique  temperature 
fields  can  be  analyzed  for  future  research. 
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THEORETICAL  MODELING  OF  THE  TRANSIENT  EFFECTS  OF  A  TOWLINE 
USING  THE  METHOD  OF  CHARACTERISTICS 


I:  Introduction 


1.1.  Background  and  Problem  Setup 

As  military  aircraft  design  has  advanced  over  the  years,  so  has  the  ability  to  detect 
and  destroy  them.  The  beginning  of  military  aviation  saw  aircraft  used  primarily  for 
scouting.  At  the  outbreak  of  World  War  I,  pilots  were  flying  around  each  other  shooting 
with  handguns  or  the  occasional  mounted  gun.  There  were  antiaircraft  defenses,  but  they 
often  missed  their  targets.  The  aircraft  themselves  were  only  mildly  effective  warfighting 
tools,  and  were  thus  not  the  focus  of  much  of  military  strategy. 

Over  the  years,  however,  military  aviation  has  become  a  vital  part  of  combat 
strategy.  The  turning  point  of  the  Pacific  Theater  of  World  War  II  was  the  Battle  of 
Midway.  A  naval  conflict,  this  battle  was  fought  purely  with  aircraft,  proving  the 
effectiveness  and  necessity  for  air  superiority.  Indeed,  at  the  end  of  the  same  war,  Japan 
found  itself  minus  two  cities  due  to  attacks  from  the  air.  One  might  wonder  what  would 
have  happened  had  Japan’s  air  defense  been  stronger  such  that  it  could  have  shot  down 
the  bombers. 

The  intense  threat  of  air  attacks  has  propelled  an  equally  intense  development  of 
antiaircraft  defenses.  From  the  introduction  of  surface  to  air  missiles  in  Vietnam  to  the 
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development  of  long  range  air  to  air  missiles  such  as  the  Phoenix,  military  aircraft  must 
increasingly  be  wary  of  threats  from  both  the  ground  and  air.  Modem  design  of  aircraft 
must  keep  these  issues  central  to  their  design. 

This  advancement  of  aircraft  design  has  meant  that  as  technology  progresses, 
much  more  care  and  work  needs  to  be  put  into  protecting  aircraft  due  to  the  immense  cost 
of  each  vehicle.  Ball  quotes  an  anonymous  source  in  his  book  as  saying,  “the  cost  of 
modem  aircraft  weapon  systems,  coupled  with  the  requirement  that  the  system  be 
effective,  makes  imperative  the  consideration  of  the  aircraft’s  survivability  throughout  the 
life  cycle  of  the  system”  (Ball,  2003:  xxvii). 

The  development  of  aircraft  for  modern  conflict  must  take  careful  consideration 
to  the  threats  that  are  posed.  This  is  the  back  and  forth  process  that  the  history  of 
technology  in  warfare  has  followed.  Antiaircraft  artillery  was  able  to  be  more  effective  at 
the  altitude  the  aircraft  were  flying,  so  aircraft  were  designed  to  fly  higher.  Heat  seeking 
missiles  were  put  into  battle,  so  chaff  and  flare  became  standard  defenses.  Radar  could 
detect  threats  and  track  aircraft,  so  stealth  aircraft  and  radar  jammers  were  designed. 

A  lot  of  the  more  modern  advances  have  been  in  electro-magnetic  warfare. 
Anything  that  uses  the  EM  spectrum  falls  under  the  scope  of  EM  warfare.  This  includes 
detection  devices  such  as  radar,  as  well  as  any  type  of  communication  device.  Control  of 
the  electro-magnetic  spectmm  in  modern  conflicts  has  become  one  the  most  vital 
guarantors  of  effectiveness. 

One  way  to  protect  against  aircraft  in  the  realm  of  EM  warfare  is  the  use  of  towed 
decoys.  A  towed  decoy  is  nothing  more  than  a  relatively  small  body  towed  behind  an 
aircraft  going  into  combat  that  can  emit  EM  signals.  It  is  connected  to  the  aircraft  by  a 
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towline  that  can  send  signals  back  and  forth,  thus  eliminating  any  need  for 
communication  through  the  air.  The  towline  also  eliminates  the  need  for  a  propulsion 
and  guidance  system  on  the  towed  body.  The  body  contains  equipment  that  can  be  used 
to  confuse  enemy  radar  into  seeing  multiple  or  different  types  of  aircraft.  This  could 
deter  an  attacker  from  attempting  an  attack,  or  force  a  threat  to  reveal  their  location  (by 
firing  a  SAM,  for  instance).  If  a  weapon  is  fired,  the  weapon  will  track  the  towed  body 
rather  than  the  aircraft.  Either  way,  the  aircraft  and  her  crew  are  protected  and  are 
allowed  to  complete  their  mission. 

One  of  the  more  common  towed  decoy  systems,  Raytheon’s  AN/ALE-50  Towed 
Decoy  System  is  currently  used  on  many  US  military  aircraft,  including  the  F/A-18E/F 
Superhornet,  the  F-16  Falcon,  and  the  B-1B  Lancer.  This  system  has  been  used  for  the 
last  decade  or  so.  More  developments  on  towed  decoys  are  underway.  A  drawing  of  this 
system  is  shown  in  Figure  1.1-1. 
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(used  without  permission)  (The  Raytheon  Company,  2006) 


Figure  1.1-1  Raytheon’s  AN/ALE-50  Towed  Decoy  System 
The  towed  body  has  the  potential  to  undergo  movement  in  combat.  It  is  often 
unwise  and  tactically  foolish  to  fly  straight  and  level  when  under  fire.  Thus,  the  towline 
and  towed  body  may  need  to  undergo  significant  maneuvers.  This  can  put  a  lot  of 
tension  on  the  towline,  and  often  causes  it  to  drift  into  the  engine  plume.  At  such  a  close 
distance,  the  towline  has  the  potential  to  be  damaged  or  severed,  thus  making  the  towed 
body  worthless. 

This  paper  seeks  to  provide  a  way  to  model  the  behavior  of  the  towed  body  and 
towline  under  transient  conditions  so  that  future  work  can  use  these  results  to  better 
design  the  system.  This  paper  also  analyzes  the  convective  transfer  of  heat  to  the 
towline,  and  develops  a  manner  in  which  to  transfer  data  through  time  and  down  the  line 
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by  using  the  method  of  characteristics.  Computer  coding  done  in  MATLAB®  is  included 
that  can  be  easily  modified  to  account  for  temperature  effects  and  transient  maneuvers. 

1.2.  Assumptions 

The  assumptions  in  this  report  were  such  that  the  analysis  of  the  system  of  towline 
and  towed  body  can  be  applied  to  almost  any  towed  body  in  any  medium.  Much  of  the 
work  in  this  paper  was  based  on  the  development  of  equations  to  model  towed  bodies 
underwater.  The  first  set  of  assumptions  involved  the  towline  behavior.  The  second  set 
involved  the  attachment  of  the  towline  to  the  body  and  the  aircraft.  The  third  set  of 
assumptions  involved  the  aerodynamic  effects.  The  fourth  set  of  assumptions  involved 
heat  transfer  effects. 

Regarding  the  towline  behavior,  the  first  assumption  was  that  compression  in  the 
line  is  impossible  (see  discussion  on  Crist’s  work,  Section  2.3).  Thus,  the  minimum 
tension  in  the  line  is  zero.  The  second  assumption  was  that  the  towline  was  inextensible. 
This  is  reasonable,  since  the  towed  body  is  also  assumed  to  be  relatively  small  and  the 
towline  is  assumed  to  be  made  of  inelastic  material.  The  third  assumption  is  that  the 
towline  is  assumed  to  be  straight  along  each  increment  (ds),  which  allows  for  numerical 
integration  at  each  point  along  the  line. 

Regarding  the  attachment  of  the  towline  to  the  body  and  the  aircraft,  the  first 
assumption  is  that  the  towline  at  its  attachment  to  the  body  and  the  aircraft  has  the  same 
behavior  as  the  body  and  the  aircraft  respectively.  This  allows  for  the  computation  of  the 
upper  and  lower  boundaries  of  the  towline.  The  second  assumption  is  that  the  towed 
body  has  no  rotational  forces  on  it.  In  other  words,  any  moments  acting  upon  it  become 
negligible.  If  the  towed  body  is  loosely  attached  (i.e.,  allowed  to  move  freely)  at  its 
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lower  end,  it  should  generally  be  parallel  to  the  freestream  velocity.  Since  many  of  these 
bodies  are  also  designed  to  have  a  streamline  shape  with  fins  (to  keep  them  in  the 
freestream),  this  assumption  is  reasonable.  The  body  is  thus  modeled  as  always  being 
parallel  to  the  freestream. 

Regarding  the  aerodynamic  effects,  the  first  assumption  is  that  the  aircraft  is 
flying  at  subsonic  speeds.  This  prevents  the  need  for  any  type  of  shock  wave  analysis. 
The  second  assumption  is  that  turbulent  effects  are  negligible,  or  rather  that  the  length 
scale  of  turbulence  is  larger  than  the  control  volume.  Thus,  the  wake  of  the  aircraft  and 
any  disturbance  from  engine  plume  do  not  come  into  play,  and  all  calculations  are  based 
on  the  assumption  of  local  laminar  flow.  The  third  assumption  is  that  skin  friction 
parallel  to  the  towline  is  negligible.  Due  to  the  problem  of  an  increasing  boundary  layer 
along  the  towline,  the  solution  for  analyzing  the  boundary  layer  becomes  tricky. 

However,  this  assumption  is  the  same  used  in  Richardson  (Richardson,  2005:20),  and 
should  not  significantly  affect  the  results. 

Regarding  the  heat  transfer  effects,  the  first  assumption  was  that  the  temperature 
of  the  medium  is  known  at  every  point  in  space  and  time.  This  can  be  provided  by 
analyzing  the  behavior  of  the  engine  plume,  the  primary  source  of  heat  disturbances.  The 
second  assumption  is  that  the  towed  body  is  sufficiently  far  behind  the  aircraft  such  that 
its  temperature  remains  constant  over  time.  This  prevents  the  need  for  extra  processes 
that  would  slow  down  the  calculations.  Since  the  primary  concern  regarding  temperature 
analysis  is  the  possibility  of  damaging  or  severing  the  towline,  there  is  no  need  to  find  the 
heat  transfer  to  the  body.  The  entire  towline  is  analyzed,  however,  all  the  way  to  the 
point  of  attachment  to  the  towed  body.  The  third  assumption  is  that  the  temperature 
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change  along  the  towline  is  gradual.  Thus,  all  heat  transfer  comes  from  the  convection 
due  to  airflow  over  the  towline,  and  any  conduction  within  the  towline  is  negligible.  The 
main  concern  is  the  maximum  temperature  in  the  line,  so  this  is  also  a  reasonable 
assumption.  The  fourth  assumption  is  that  air  pressure  is  constant.  Thus,  values  of  air 
density,  viscosity,  and  thennal  conductivity  are  purely  functions  of  air  temperature. 

1.3.  Preview  of  Paper 

There  are  5  chapters  in  this  thesis,  the  first  one  (here)  giving  an  introduction  to  the 
project’s  topic.  The  rest  of  the  paper  is  as  follows: 

Chapter  2  contains  a  summary  of  background  reading  that  was  utilized  in 
researching  this  topic.  Since  this  paper  is  not  experimental,  the  previous  work  became 
vital  in  understanding  the  development  of  techniques  to  model  towline  behavior.  The 
chapter  concludes  with  a  description  of  the  method  of  characteristics,  a  vital  aspect  of  the 
methodology. 

The  methodology  is  contained  in  chapter  3.  The  development  of  the  governing 
equations,  as  well  as  converting  the  equations  of  motion  using  the  method  of 
characteristics  is  explained  in  great  detail.  A  set  of  equations  is  determined  that  can 
model  the  behavior  of  the  towline  under  different  perturbations  with  given  initial  data.  In 
addition,  a  methodology  to  analyze  heat  transfer  is  explored,  and  equations  are  developed 
to  model  heat  transfer  in  a  predetermined  temperature  field.  These  equations  can  be  used 
with  the  aid  of  computer  software,  such  as  MATLAB®  (used  for  this  paper),  to  solve  for 
the  behavior  of  a  towline  under  transient  effects. 

Results  of  the  computer  simulation  are  contained  in  chapter  4.  Some  example 
scenarios  are  offered  and  analyzed,  including  steady  state  behavior,  single  axis 
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perturbations,  a  multiaxial  perturbation,  and  heat  transfer  to  a  moving  towline  under 
constant  ambient  temperature.  These  results  were  compared  to  previous  work. 

Chapter  5  offers  general  conclusions  to  the  research.  It  also  lays  out  some 
suggestions  for  future  work,  including  the  motivation  for  such  work. 
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II:  Literature  Review 


A  lot  of  past  research  has  been  done  on  towed  lines  and  is  included  here  as  a 
source  for  background  information.  Of  particular  note  are  the  works  of  Richardson,  as 
well  as  Schram  and  Reyle.  Most  of  the  research  in  this  paper  is  based  on  their 
methodologies.  The  last  section  of  the  literature  review  will  include  a  description  of  the 
method  of  characteristics. 

2.1.  Richardson 

“Parametric  Study  of  the  Towline  Shape  of  an  Aircraft  Decoy,”  completed  2005, 
develops  a  method  to  compute  the  steady  state  form  of  a  towline  behind  and  aircraft  in 
flight.  This  thesis,  completed  by  ENS  Tyler  Richardson  (USNR),  is  the  foundational 
work  that  this  paper  is  continuing. 

Richardson  derived  equations  of  motion  by  taking  into  account  airflow  across  a 
towline  with  a  towed  body  (causing  drag),  and  the  weight  of  both  the  towline  and  the 
body.  These  equations  were  manipulated  such  that  seven  governing  differential 
equations  could  be  used  to  describe  the  behavior  of  the  towline  from  the  body  up  to  the 
aircraft.  All  equations  were  with  respect  to  the  arc  length  of  the  towline,  were 
nondimensionalized,  and  could  thus  be  integrated  along  the  towline  (see  discussion  on 
the  method  of  characteristics,  Section  2.8).  These  equations  involved  four  second  order 
ODE’s,  three  of  which  described  position  (three  axes)  and  one  of  which  described 
tension,  and  three  first  order  ODE’s  that  described  position  (three  axes).  Essentially,  only 
four  equations  were  uniquely  derived  by  Richardson,  noting  that  the  first  order  position 
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terms  become  the  integrals  of  the  second  order  terms  (i.e., 


dx 

~dl 


r  d2X 


dl ,  where  dl  is  the 


change  in  length  along  the  towline).  These  seven  governing  equations  were  placed  into 
the  ode45  solver  in  MATLAB®  and  integrated  along  the  towline  from  the  body  up  to  the 
aircraft.  This  procedure  was  repeated  for  different  parameters  (velocity,  weight,  etc.)  in 
order  to  determine  the  steady  state  behavior  of  the  line  under  various  conditions. 

The  most  important  results  from  this  work  include  MATLAB®  code  that  takes 
initial  data  and  outputs  the  position,  velocity,  and  tension  of  a  towline  in  steady  state 
flight.  Richardson’s  work  was  compared  to  real  data  and  detennined  to  be  accurate. 

Thus  his  code  will  serve  as  a  comparison  here.  Richardson’s  code  was  also  used  in  this 
paper  to  set  up  the  initial  conditions  for  transient  behavior.  In  addition,  his  work  serves 
as  a  reference  for  developing  the  method  by  which  to  solve  the  problem  at  hand. 

Although  similar  techniques  were  used  for  this  paper,  one  primary  difference  that 
Richardson’s  thesis  utilized  was  the  ode45  solver  in  MATLAB®.  The  method  for  the 
research  of  this  paper  does  not  allow  for  continuous  integration,  however,  since  a  fixed 
grid  system  is  to  be  set  up  in  advance  for  numerical  integration.  Richardson  analyzes  a 
steady  state  position  of  the  towline  before  he  analyzes  its  behavior.  This  paper  is, 
however,  concerned  with  the  transient  behavior,  and  as  a  result  deals  with  changing 
variables  over  time.  Thus,  Richardson  converts  all  of  his  equations  to  functions  of  one 
independent  variable,  whereas  this  work  requires  the  use  of  two  (see  discussion  on  the 
method  of  characteristics  in  Section  2.8).  This  addition  of  the  time  element  precludes  the 
use  of  the  ode45  solver.  Richardson  notes  himself  in  his  conclusions  and  remarks  that  “a 
numerical  integration  routine  for  the  system  of  differential  equations  needs  to  be  written 
for  the  transient  case”  (Richardson,  2005:51). 
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2.2.  Schram  and  Reyle 


“A  Three-Dimensional  Dynamic  Analysis  of  a  Towed  System,”  published  in 
1968,  explores  the  behavior  of  an  underwater  towed  body  behind  a  ship.  Using  a  purely 
theoretical  model,  the  paper  starts  off  by  deriving  the  equations  of  motion  from  Newton’s 
second  law  in  three  dimensions.  These  take  into  account  not  only  tension  from  the  ship, 
but  include  hydrodynamic  towline  forces.  A  coordinate  transformation  was  applied  so 
that  all  values  were  respective  of  the  towline  itself.  For  boundary  conditions,  they 
assumed  that  “the  towline  at  its  point  of  attachment  to  the  towing  ship  must  have  the 
same  motion  as  the  ship  . . .  [and]  the  towline  must  have  the  same  motion  as  the  body  at 
its  point  of  attachment”  (Schram  and  Reyle,  1968:216).  They  also  set  the  tow  point  to 
“coincide  with  the  center  of  mass  of  the  body”  in  order  to  eliminate  pitch,  roll,  and  yaw 
effects  (Schram  and  Reyle,  1968:216).  In  order  to  understand  the  dynamic  effects  of 
towline  motion,  the  method  of  characteristics  was  applied  such  that  all  values  vary  solely 
with  length  and  time.  They  finished  by  developing  a  numerical  procedure  to  solve  the 
dynamic  solution. 

Schram  and  Reyle  made  a  few  conclusions  that  are  important  for  this  paper.  First, 
they  found  that  “the  speed  at  which  transverse  disturbances  are  propagated  in  the  x-y  and 
the  z-y  planes  are  the  same”  (Schram  and  Reyle,  1968:217).  Second,  they  found  that  “if  a 
towline  is  not  straight,  each  type  of  disturbance,  transverse  or  longitudinal,  influences  the 
other”  (Schram  and  Reyle,  1968:219).  In  other  words,  a  solely  longitudinal  disturbance 
could  cause  transverse  motion  farther  down  the  towline  as  well  and  vice  versa.  Lastly, 
they  found  that  the  transfer  function,  which  is  “the  ratio  of  the  resultant  amplitude  of  the 
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body  motion  to  the  amplitude  of  the  motion  of  the  ship”  (Schram  and  Reyle,  1968:220) 
decreases  if  either  the  towline  length  or  towing  speed  increases.  With  a  longer  length, 
disturbances  have  a  longer  distance  to  travel,  thus  being  subject  to  a  greater  amount  of 
damping.  With  a  faster  speed,  the  towline’s  pitch  angle  is  reduced,  causing  a  greater 
percentage  of  disturbances  to  be  transversal.  Transverse  disturbances  were  found  to  be 
damped  significantly  more  than  longitudinal  disturbances. 

The  research  was  summarized  in  a  paper  published  in  the  Journal  of  Hydronautics 
in  1968,  which  was  used  as  a  primary  point  of  reference  here.  This  paper  was  essentially 
a  summarization  of  a  PhD.  thesis  completed  earlier  the  same  year  (Reyle  was  the  PhD. 
advisor)  at  Rutgers  University.  Both  the  journal’s  paper  and  the  PhD.  thesis  were  used  in 
conjunction  to  understand  their  methodology. 

Schram  and  Reyle’s  research  provides  the  basis  of  much  of  the  methodology  in 
this  paper.  Their  method  will  not  be  explained  in  more  detail  here  due  to  the  complexity 
of  it.  However,  the  methodology  section  of  this  paper  explains  much  of  what  they  did. 
Their  application  of  the  method  of  characteristics  in  modeling  the  behavior  of  the  towed 
body  has  been  relied  upon  heavily.  In  addition,  their  conclusions  for  towline  behavior 
will  serve  as  a  comparison.  The  results  of  this  research  will  differ  in  the  sense  that 
Schram  and  Reyle  studied  underwater  behavior  behind  a  towing  body  at  constant  altitude 
and  ignored  heat  transfer  and  temperature  effects.  This  paper  will  explore  the  behavior  of 
the  towed  body  in  air  at  differing  altitudes  and  set  up  a  method  to  analyze  heat  transfer 
and  temperature  effects. 
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2.3.  Crist 


“Analysis  of  the  Motion  of  a  Long  Wire  Towed  from  an  Orbiting  Aircraft,” 
published  in  1970,  studies  the  behavior  of  a  long  cable  towed  behind  an  aircraft 
undergoing  oscillations,  to  include  wind  shear  effects  on  cable  motion.  The  paper  relies 
upon  previous  steady  state  research  on  the  T  AC  AMO  (TAke  Charge  And  Move  Out) 
cable.  It  assumes  a  lumped  mass  model  of  the  cable  and  derives  Lagrange’s  equations  of 
motion.  The  equations  were  solved  numerically  for  vertical  oscillations  of  a  towing 
aircraft  in  constant  radius  and  altitude,  as  well  as  transition  from  orbit  to  straight  and 
level  flight.  Crist  gives  the  method  for  wind  shear  analysis,  but  does  not  provide  results. 

Crist’s  concludes  that  “the  effects  of  slack  cable  must  be  included  in  most 
analyses,”  that  “the  tension  at  the  drogue  can  be  considerably  higher  than  the  equilibrium 
tension  when  the  aircraft  oscillates  vertically,”  and  that  the  “tension  in  the  cable  during 
transition  from  orbit  can  be  considerable,  even  for  a  smooth  transition  from  orbit”  (Crist, 
1970:73).  Perhaps  more  importantly,  he  notes  that  “compressive  force  in  the  cable 
segment...  is,  of  course,  impossible  since  everyone  knows  ‘you  can’t  push  on  a  rope’” 
(Crist,  1970:65).  He  also  developed  an  analysis  to  calculate  the  motion  of  a  towed  cable 
behind  an  orbiting  aircraft  following  a  prescribed  flight  path.  Due  to  the  high  speed  of 
aircraft,  however,  slack  conditions  will  not  be  analyzed  in  this  paper,  since  they  create 
negative  tension  errors.  More  on  this  will  be  discussed  in  the  conclusions. 

Crist’s  analysis  of  a  towing  cable  undergoing  aircraft  motion  will  be  used  as  a 
source  for  comparison  for  this  paper.  Due  to  computing  power  at  the  time  his  paper  was 
written,  a  discrete  mass  approach  was  used  in  the  analysis.  Crist  notes,  however,  that  a 
good  approximation  can  be  made  to  a  continuous  mass  approach  with  only  a  few  mass 
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points  (Crist,  1970:61).  Since  this  paper  is  using  numerical  integration  along  a  towline,  it 
will  also  use  the  discrete  mass  approach,  although  more  accurate  results  will  be  available 
due  to  much  faster  computing  time.  In  addition,  this  paper  will  develop  methods  to 
analyze  a  whole  range  of  motion,  not  just  vertical  oscillations  or  transition  to  straight  and 
level  flight.  As  in  the  Schram  and  Reyle  paper,  Crist  does  not  consider  heat  transfer  or 
temperature  effects  either. 

2.4.  Chapman 

“Towed  Cable  Behaviour  During  Ship  Turning  Manoeuvres,”  published  in  1984 
explores  the  behavior  of  an  underwater  towed  body  in  large  and  small  radius  turns.  In 
large  radius  turns,  the  system  behaves  similarly  to  a  steady  state  straight-tow  behavior. 
However,  in  small  radius  turns,  the  system  collapses,  causing  the  towed  body  to  drop  in 
depth  (no  cable  tension),  followed  by  a  quick  tug  (high  cable  tension).  This  paper’s 
purpose  was  to  detennine  a  non-dimensionalized  minimum  radius  at  which  collapse  can 
be  avoided. 

Chapman  concluded  that  “in  a  gradual  turn  (with  a  relatively  slow  rate  of  turn)  the 
cable  moves  towards,  and  usually  attains,  an  equilibrium  configuration  that  is  a  slightly 
perturbed  fonn  of  the  profile  adopted  during  straight  towing.  The  [towed  body]  moves 
towards  a  circular  course  which  is  concentric  with  the  ship’s  course  but  of  slightly 
smaller  radius.  As  a  consequence  the  speed  of  the  [towed  body]  is  reduced  resulting  in  a 
slight  increase  of  depth”  (Chapman,  1984:  353).  Also,  he  concludes  that  in  a  sharp  turn, 
the  straight-tow  profile  collapses,  causing  the  towed  body  to  drop  sharply.  This  results  in 
increased  tension  at  the  top  of  the  cable,  and  could  cause  slackness  near  the  towed  body. 
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Thus,  the  system  will  most  likely  never  reach  equilibrium  in  the  turn.  He  found  that  the 
transition  from  equilibrium  turns  to  a  collapsed  system  happens  over  a  very  small  range 
of  radii,  thus  giving  a  fairly  defined  minimum  turning  radius. 

Chapman’s  analysis  and  methods  will  be  used  in  this  paper  to  help  understand 
small  radius  turns.  His  analyses  involved  an  underwater  body,  which  will  exhibit 
significantly  more  drag  than  a  body  in  air.  Also,  ship  speeds  will  be  significantly  slower 
than  those  of  aircraft.  Thus,  while  the  concern  of  a  collapsed  system  in  sharp  turns  will 
be  considered,  it  will  most  likely  not  be  an  issue  except  in  extreme  cases.  Due  to  the 
restriction  in  this  paper  of  no  slack,  Chapman’s  research  is  referenced  here  since  it  can  be 
used  to  prevent  slack  conditions  for  future  work. 

2.5.  Dowling 

“The  Dynamics  of  Towed  Flexible  Cylinders  Part  I:  Neutrally  Buoyant 
Elements,”  published  in  1987,  investigates  the  transverse  vibrations  of  a  thin,  flexible 
neutrally  buoyant  cylinder  under  nominally  constant  towing  conditions.  The  cylinder  is 
assumed  to  have  a  very  small  bending  stiffness,  and  stretches  out  due  to  viscous  forces 
along  its  length. 

Dowling  focuses  primarily  on  a  flexible  towed  body,  whereas  this  paper’s 
research  is  concerned  with  inflexible  towed  bodies.  Also,  he  investigates  underwater 
behavior  with  a  towing  body  of  constant  altitude.  This  paper  will  be  used  primarily  to 
develop  methods  for  forces  on  the  towed  body. 
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2.6.  Ames 


Nonlinear  Partial  Differential  Equations  in  Engineering,  published  in  1965,  gives 
a  detailed  understanding  of  the  method  of  characteristics,  and  provides  a  foundation  for 
the  mathematical  understanding  behind  the  derivations  used  in  this  paper.  Specifically, 
this  book  helps  set  up  the  mesh  scheme  and  the  equations  to  carry  data  from  one  timestep 
to  the  next  for  the  method  of  characteristics,  which  Ames  notes  is  based  on  a  previous 
scheme  set  up  by  Courant,  Isaacson,  and  Rees.  Ames  notes  that  this  mesh  is  based  on  a 
first-order  scheme  (Ames,  1965:446),  which  becomes  significant  in  the  error  analysis  of 
this  work.  This  method  is  used  in  order  to  avoid  the  “messy  two-dimensional 
interpolation  in  the  characteristic  grid”  (Ames,  1965:444)  by  setting  up  the  mesh  points  in 
advance  and  interpolating  as  the  computation  proceeds,  making  the  interpolation  one¬ 
dimensional. 

Ames  also  proposes  a  “hybrid”  method  by  Hartree,  which  provides  a  second-order 
truncation  error.  This  method  could  be  useful  for  future  analysis,  but  is  not  used  here  due 
to  the  iteration  requirement  to  solve  the  system  of  equations  at  each  timestep.  The 
method  based  on  Courant,  Isaacson,  and  Rees  requires  no  iteration,  and  is  the  same 
method  used  in  Schram  and  Reyle.  Thus,  this  method  will  be  used  for  the  analysis  of  this 
paper. 

Schram  and  Reyle  relied  heavily  on  the  work  of  Ames,  and  this  paper  will  follow 
in  their  footsteps.  It  will  be  used  as  a  point  of  reference  for  much  of  the  math,  and  will 
help  with  the  derivation  of  any  new  equations. 
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2.7.  Tannehill,  Anderson,  and  Pletcher 


Computational  Fluid  Mechanics  and  Heat  Transfer,  published  in  1997,  serves  as  a 
foundational  basis  for  a  lot  of  CFD  work  today.  Since  this  paper’s  work  depends  on  the 
stepping  method  of  solving  a  system  of  equations,  the  CFD  book  is  a  vital  resource  in 
determining  numerical  differencing  schemes,  as  well  as  understanding  the  method  of 
characteristics.  This  book,  like  Ames,  will  be  used  primarily  as  a  reference  tool  to 
develop  the  governing  equations. 

2.8.  Method  of  Characteristics 

As  alluded  to  earlier,  both  Ames  as  well  as  Tannehill,  Anderson,  and  Pletcher 
give  detailed  descriptions  of  the  method  of  characteristics.  These,  in  conjunction  with 
Schram  and  Reyle,  were  used  to  understand  the  method  of  characteristics.  This  method 
will  be  explained  in  its  general  principles  here. 

The  method  of  characteristics  is  a  technique  that  is  used  to  solve  nonlinear 
hyperbolic  partial  differential  equations.  What  this  method  does,  in  essence,  is  convert 
the  equations  of  a  system  into  new  equations  that  are  constant  along  characteristic  lines. 
In  the  application  here,  the  system  of  equations  is  rewritten  to  be  a  function  of  two 
independent  variables,  which  can  be  held  constant  along  the  characteristic  lines. 

An  analogy  that  can  be  used  to  describe  this  is  a  wave  on  the  ocean  traveling  in 
space  and  time.  From  shore,  looking  straight  along  a  line  in  the  water,  an  observer  sees 
the  water  move  up  and  down  in  a  different  manner  at  every  point  along  that  line  (Figure 
2.8-1).  Thus,  the  behavior  is  a  function  of  time  and  location  (i.e.,  position  on  the  line  the 
observer  is  concerned  about).  If  the  reference  frame  is  changed,  however,  such  that  the 
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observer  is  on  the  top  of  one  wave,  with  his  line  of  sight  along  the  top  of  that  wave,  the 
behavior  will  look  the  same  at  all  time  and  position  along  his  line  of  sight  (Figure  2.8-2). 
The  method  of  characteristics  uses  this  type  of  technique  in  order  to  “march”  data  from 
one  timestep  to  the  next  by  finding  characteristic  lines  on  which  the  solution  is  constant. 
Applying  the  wave  analogy  to  the  towline  problem,  the  towline  at  a  certain  time  becomes 
a  waveline  in  its  direction  of  motion  (left  to  right  in  figure),  and  the  change  in  time  is 
represented  by  each  successive  waveline  (top  to  bottom  in  figure,  and  not  to  be  confused 
with  the  time  described  in  the  analogy). 


Reference  line  fixed  in  space 


Location  of  observer  -  fixed  on  shore 

Figure  2.8-1  Initial  Reference  Frame 
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Reference  line  fixed  on  wave 


\  Location  of  observer  -  fixed  on  wave 


Figure  2.8-2  New  Reference  Frame  -  fixed  to  wave 


This  method  is  vital  for  finding  a  solution  to  the  problems  of  this  paper.  What 
will  be  shown  later  is  the  development  of  six  PDE  equations  of  motion  with  respect  to 
time  and  length  along  the  towline.  In  order  to  get  a  total  solution  to  the  system,  these 
equations  will  be  converted  into  characteristic  equations,  solved  for  all  positions  along 
the  towline  for  one  timestep,  then  “marched”  to  the  next  timestep  using  the  characteristic 
lines.  To  follow  the  wave  analogy  represented  in  Figure  2.8-2,  the  information  at  one 
point  on  the  towline  is  carried  to  the  next  timestep  at  a  different  position  on  the  towline. 
This  new  position  is  detennined  by  the  characteristic  line  and  the  change  in  time. 
Another  way  to  show  this  is  in  Figure  2.8-3,  which  shows  a  solution  in  time  and  space 
that  will  be  used  for  this  paper,  where  s  is  the  length  along  a  towline  and  t  is  time.  This 
figure  is  based  on  one  devised  by  Tannehill,  Anderson,  and  Pletcher  (Tannehill, 
Anderson,  and  Pletcher,  1997:361),  but  could  be  found  in  most  detailed  descriptions  of 
the  method  of  characteristics. 
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- Characteristic  Lines 

•  Solution 

Figure  2.8-3  Characteristic  Net 

Note  that  the  solution  can  be  found  at  every  intersection  of  the  characteristic  lines. 
Thus,  new  data  is  available  at  every  interval.  Since  the  characteristic  lines  carry  constant 
data  along  their  entire  path,  middle  steps  can  be  skipped  if  desired.  If  new  lines  were 
started  from  every  5  and  t  position,  the  density  of  the  mesh  would  double.  In  fact,  the 
mesh  size  can  be  reduced  to  be  as  small  as  desired  by  noting  that  any  new  used  data  point 
can  transfer  information  along  newly  created  characteristic  lines  from  that  point. 

The  characteristic  net  is  defined  by  the  characteristic  lines.  The  positive  sloped 
characteristic  lines  all  have  the  same  slope,  but  start  at  equally  spaced  positions  along  the 
5-axis.  Similarly,  the  negative  sloped  lines  all  have  the  same  slope,  starting  at  the  same 
positions  as  their  counterparts.  Due  to  the  hyperbolic  nature  of  the  equations  used  in  the 
method  of  characteristics,  these  two  slopes  are  different  only  in  sign,  and  are  found  by 
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finding  the  roots  of  the  system  of  equations.  These  roots  are  called  the  “characteristic 
lines,”  and  make  up  the  “characteristic  net”  shown  in  Figure  2.8-3. 

Thus,  the  application  of  the  method  of  characteristics  necessitates  first  that 
characteristic  equations  are  derived  such  that  they  remain  constant  along  the 
characteristic  lines.  This  is  done  by  converting  a  set  of  hyperbolic  PDE’s  into  a  set  of 
ordinary  differential  equations  where  each  ODE  comes  from  two  separate  PDE’s.  Each 
new  ODE  is  then  integrated  along  the  characteristic  lines.  Ames  (Ames,  1965:  422)  notes 
that  these  calculations  must,  in  most  cases,  be  carried  out  numerically.  Exact  integration 
by  this  method  is  very  rare.  This  paper  will  use  the  method  of  characteristics  with 
numerical  integration. 
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Ill:  Methodology 


The  subsequent  derivations  for  towline  motion  come  from  an  adaptation  of  those 
in  Schram  and  Reyle.  While  their  work  is  heavily  relied  upon  here,  they  were  not  always 
followed  directly.  Effort  has  been  made  to  note  when  their  work  was  followed.  The  heat 
transfer  derivations  were  completely  independent  of  any  past  work. 

3.1.  Derivation  of  Governing  Equations 

The  governing  equations  and  their  derivations  rely  very  heavily  upon  work  done 
by  Schram  and  Reyle  (1968).  Similar  notation  and  procedure  will  be  used,  and  a  great 
deal  of  effort  has  been  used  to  reduce  any  ambiguity  in  their  work,  and  to  make  the 
procedure  used  in  this  paper  as  clear  as  possible. 

Axis  Transformation. 

In  order  to  derive  a  system  of  equations  for  use  with  the  method  of  characteristics, 
all  variables  will  need  a  manner  to  convert  between  coordinate  axes  aligned  with  space 
(aircraft  at  steady  state),  and  those  aligned  with  the  towline.  Thus,  a  matrix  [A]  must  be 
found,  which  will  allow  for  the  transformation  of  components  (position,  velocity,  etc.) 
between  space  and  towline  axes.  This  transformation  is  shown  in  Figure  3.1-1,  Figure 
3.1-2,  and  Figure  3.1-3 
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Figure  3.1-1  Spatial  orientation  of  x-y  plane 


Figure  3.1-2  Spatial  orientation  of  z-y  plane 


Figure  3.1-3  Spatial  orientation  of  z-x  plane 

where  (j)  is  the  azimuth  angle  and  6  is  the  polar  angle.  A  note  should  be  made  about  6  in 
regards  to  its  appearance  in  the  diagrams.  This  angle  refers  to  a  second  transform  about 
the  towline  X-axis,  and  the  two-dimensional  silhouette  of  Figure  3.1-2  and  Figure  3.1-3 
do  not  show  the  fact  that  9  can  go  in  or  out  of  the  page  as  well.  A  better  visual 
representation  is  in  Figure  3.1-6,  while  making  sure  to  note  the  transitional y  ’-axis. 

The  spatial  axes  are  oriented  with  the  v-axis  aligned  directly  vertical,  the  x-axis 
aligned  with  the  direction  of  the  aircraft,  and  the  z-axis  aligned  perpendicular  to  the 
aircraft’s  line  of  flight.  The  towline  axes  are  oriented  with  the  7-axis  aligned  directly 
along  the  towline,  the  X-axis  aligned  downward  from  the  towline  and  in  the  x-y  plane, 
and  the  Z-axis  aligned  perpendicular  to  the  towline  and  the  X-axis.  Figure  3.1-4  sets  up  a 
three  dimensional  plot  for  transformation. 
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The  first  transformation  changes  the  spatial  coordinates  into  the  towline  X-axis 
and  a  mid-step  y  ’-axis  by  rotating  about  the  z-axis.  This  is  shown  in  Figure  3.1-5.  The 
angle  (f)  is  used  here  and  exists  only  in  the  spatial  x-y  plane.  Thus,  (f>  is  not  hard  to 
represent  in  the  other  diagrams  of  this  paper. 


Figure  3.1-5  First  Transformation 
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Note  that  <f>  is  defined  as  positive  up  from  the  x-axis  such  that  a  zero  value  lines  up  they  ’ 
andx  axes. 

From  here,  the  following  are  derived: 

/(^O  =  f(x)  sin  </>  ~  f(y)  cos  (/)  1.1 

/  (y ')  =  /' (*)  cos  </>  +  f  ( v)  sin  (j)  1.2 

where  f{  )  refers  to  any  variable  in  each  respective  axis.  Note  that  the  y  ’-axis  has  no  real 
value  except  for  its  use  in  the  process  of  coordinate  transfonnation.  As  will  be  seen  in 
the  next  transformation,  however,  a  6  angle  of  zero  makes  the  y  ’  and  y  axes  line  up,  and 
Figure  3.1-5  is  sufficient  to  model  the  axis  transformation. 

At  this  point,  the  transfonnation  from  space  to  the  X  towline  component  is 
complete.  The  second  transformation  will  provide  a  solution  for  the  Y  and  Z  components 
in  the  towline  by  rotating  about  the  new  X-axis.  This  transformation  is  shown  in  Figure 
3.1-6.  As  an  aside,  the  angle  6,  used  for  the  second  transformation,  exists  in  the  z-y  ’ 
plane,  where  y  ’  is  the  transitory  axis.  Thus,  it  is  hard  to  represent  accurately  in  many  of 
the  diagrams. 


26 


z 


$ 

Z‘ 

Figure  3.1-6  Second  Transformation 
Thus,  the  following  are  derived: 

/(7)  =  /  (y ')  cos  B  +  /(z)  sin  e 

=  f(x)cos<ficos8  +  /( v)  sin  <j>  cos  6  +  /(z)sin# 

/(Z)  =  -/( v  ’)  sin  9  +  /(z)  cos  0 

=  -/ (x)  cos  ^  sin  0-  f  ( v)  sin  (f>  sin  6  +  f  (z)  cos  61 

From  these  equations,  a  matrix  [A]  is  generated  such  that  all  axis  transformations 

can  be  computed  in  a  simple  manner.  By  taking  the  components  of  x,  y,  and  z,  one 

generates  the  matrix  as  follows: 


[A]  = 


sin^ 

cos^cos# 

-cos^sin# 


-cos  <9  0 

sin  (f)  cos  9  sin  9 

-  sin  (j)  sin  0  cos  0 


1.5 


This  matrix  can  be  used  to  transform  from  space  to  towline,  and  vice  versa.  This  is 
accomplished  quite  easily,  and  is  printed  here  for  reference  (to  prevent  confusion,  TL  and 
S  will  denote  towline  and  space  respectively): 
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1.6 


~f{X) 

fix) 

f(Y ) 

=  [A\ 

f(y ) 

J(Z)\ 

TL 

_/(*)_ 

[/(*)  f(y)  /(*)]*  =[/W  f(X)  fV)]a[A  L7 

Now  that  a  manner  in  which  to  transfonn  between  the  spatial  and  towline  axes  has  been 
determined,  the  governing  equations  may  be  developed. 


Equations  of  Motion. 

Since  this  system  contains  six  unknowns  (U,V,W,  <f> ,  0,  T)  in  the  towline  reference 
frame  at  any  point  in  time  and  space,  six  equations  of  motion  will  be  needed  to  describe 
the  entire  behavior  of  the  towline. 

Using  Newton’s  Second  Law  of  Motion  (flipped  from  standard  notional  for 
convenience) 

ma  =  F  1.8 


and  applying  it  to  the  free  body  diagram  of  the  towline  in  its  spatial  coordinates  (Figure 
3.1-4),  Schram  cites  Cristescu  (Schram,  1968:6)  with  determining  that  the  equations  of 
motion  for  an  elastic  line  become: 


d2x 


dt  ds 


T  dx 
(1  +  e)  ds 


+  F 


d2y  _  d 
dt 2  ds 


l  +  F  -Wt 
(1  +  e)  ds  ' 


d2z  d 


dt  ds 


T  dz 
(1  +  e)  ds 


+  F 


It  should  be  noted  that  these  tenns  are  all  per  meter  of  towline  length. 


1.9 

1.10 


1.11 
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In  these  equations,  tension  is  divided  by  a  percent  change  in  extension  of  the 
towline.  Assuming  an  inextensible  towline,  which  this  paper  will  assume  just  as  Schram 
and  Reyle  did,  the  extension  terms  drop  out.  The  8x18  s,  8yl  8  s,  and  8  zl  8  s  terms  are 
used  to  project  tension  into  the  respective  axes,  since  tension  is  purely  directed  along  the 
towline  (the  5-direction  is  used  here,  which  is  interchangeable  with  the  7- axis).  Through 
a  simple  application  of  the  matrix  transformation  represented  in  Equation  1 .7,  these  tenns 
become: 


dx 
8 s 

8y_ 

8s 


-cos#  cos  ^ 


-cos#  sin  <f> 


8z 

8s 


=  -sin# 


where  the  negative  term  is  due  to  the  orientation  of  s  being  positive  down  the  towline 
towards  the  body. 

Substituting  Equations  1.12  through  1.14  into  Equations  1.9  through  1.11, 
replacing  the  second  derivative  (acceleration)  terms  with  equivalent  velocity  terms,  and 
rearranging  slightly,  these  equations  reduce  to: 


8u  _  1 
8t  // 


—(-T  cos#cos^)  +  F 

8s y  ’  x 


8v  _  1 
8t  fd 


—  (-T  cos#sin^)  +  Fv  -  Wt 
8s 


8w 

8t 


1 

A 


8_ 

8s 


[-T  sin#)  +  F. 


1.12 

1.13 

1.14 


1.15 

1.16 

1.17 
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Taking  the  three  towline  acceleration  terms  and  applying  them  to  the  axis 
transformation  in  Equation  1 .6  yields  the  equations  of  motion  in  towline  form,  noting 
first  the  velocity  transformations: 


U  =  u  sin  (f)  -  v  cos  <f) 

V  =  i/cos^cos6,  +  vsin0cos#  +  wsin# 


W  =  -u  sin#cos^-vsin6,sin0  +  wcos9 


Thus,  the  equations  of  motion  become: 


dU 

dt 


du  .  dv 

—  sin® - cos  ® 

dt  dt 


+  (ncos^  +  vsin^)  — 
V  ’  dt 


dV_ 

dt 


f  du  ,  dv  . 

—  cos  ®  H - sin® 

dt  dt 


A  dw 

cosd  +  —  sin^ 
j  dt 


dd)  dd 

+  ( v  cos  <j)  -  u  sin  0)  cos  9  — —  ({u  cos  <j)  +  v  sin  </>)  sin  9  +  wcos9)  — 


dW 

dt 


du 


dv  . 


■  n  dw  „ 
sin6^  +  —  cos  9 
dt 


—  cos^H - sinizi 

dt  dt  j 

■  ( -u  sin  <f)  +  v  cos  (f> )  sin  9  -  (( w  cos  (f>  +  v  sin  cos  9  +  w  sin  9 ) 


The  bracketed  terms  come  from  the  direct  transfonnation  of  the  U,  V,  and  W 
towline  variables  (Equations  1.18,  1.19,  and  1.20),  while  the  extra  terms  come  from  the 
derivative  due  to  the  product  rule.  By  utilizing  Equation  1.7,  the  spatial  velocity  terms 
can  also  be  transformed  into  their  towline  forms: 

u  =  U  sin  (f)  +  V  cos  (f>  cos  9  -  W  cos  (f>  sin  9 
v  =  —U  cos  (j)  +  V  cos  9  sin  <f)  -  W  sin  9  sin  (j) 
w  =  Vsin9  +  W  cos  9 


1.18 

1.19 

1.20 

1.21 

1.22 


1.23 


1.24 

1.25 

1.26 
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The  spatial  forces,  Fx,  Fy,  and  Fz,  are  also  transformed  in  the  same  manner  to  the 
towline  forces  Fx,  Fy,  and  Fz.  Combining  the  above  equations  by  replacing  all  of  the 
spatial  terms  with  towline  terms,  everything  can  now  be  transformed  into  towline 
parameters.  Thus,  the  equations  of  motion  in  towline  axes  become: 

First  Three  Equations  of  Motion: 


M 


8U 

dt 


+  {V  cos  0  -W  sin  6) 


d(j) 

dt 


-T  cos  6 


d(j) 

ds 


Fx  -  Wt  cos  (f) 
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F 


dV_ 

dt 


+  U  cos  0 


d(j) 

dt 


- vFy  -  W,  sin ^ cos# 

ds 


1.28 


F 


d(j)  d6  dW 
U  sin  9  —  -  V - 


dt 


dt  dt 


dO 

=  T - Fz-Wt  sin  (j)  sin  6 

ds 


These  equations  will  be  useful  later  in  the  method  of  characteristics,  and  provide 
the  first  3  equations  of  motion.  All  partial  derivatives  are  with  respect  to  either  time  or 
arc  length  ( s ),  which  is  valuable  since  the  two  independent  variables  used  in  the  method 
of  characteristics  will  be  time  and  arc  length. 

The  next  three  equations  can  be  found  by  simply  taking  the  partial  derivative  of 
the  towline  accelerations  with  respect  to  the  arc  length  of  the  towline.  These  equations 
become  exactly  the  same  as  Equations  1.21,  1.22,  and  1.23,  except  that  the  partial 
derivatives  are  now  in  terms  of  s  instead  of  t: 


dU 

ds 


du  .  dv 
—  sin® - cos® 

ds  ds 


d(j)  .  d(f> 

+  u  cos  (p - 1-  v  sin  (h  — 

ds  ds 


1.29 


1.30 


dV_ 

ds 


du 


dv 


A 


—  cos^H - sin<zi 

ds  ds 


n  dw  . 
cos#  +  — sin  0 

ds 


+  ( v  cos  (j)  -  u  sin  (j) )  cos  6  —  -  ( (u  cos  (f>  +  v  sin  <ft)  sin  6  +  w  cos  0 )  — 

ds  ds 


1.31 
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8W 

ds 


du 


dv 


—  cos  (j)  H - sin^ 

yds  ds  j 


■  n  dw  n 

sin#  + —  cos  # 

ds 

W</> 


-  ( -u  sin  </>  +  v  cos  (j) )  sin  6  ~[{u  cos  </>  +  v  sin  (j>)  cos  #  +  w  sin  6 ) 

The  partial  derivatives  of  the  spatial  velocities  with  respect  to  arc  length  can  be  easily 
related  in  the  following  manner  with  the  use  of  Equations  1.12,  1.13,  and  1.14: 


(  Piv-A 


ds  ds 


yvi  j 


yvs  j 


=  —  (-cos#  COS  0) 
dt 


=  — (-cos#  sin  (f>) 
dt 


dw  d 


f  dz^ 


ds  ds 


ydtj 


f  dz^ 


dt 


\ds; 


—  (-sin#) 
dt 


By  combining  these  with  the  velocity  relationships  from  Equations  1.24,  1.25,  and  1.26, 
one  derives  the  final  three  equations  of  motion  in  towline  form: 

Last  Three  Equations  of  Motion: 

=  —  cos  #  +  (V  cos  #  -  W  sin  #)  — 
ds  dt  ds 


dV  dO  d(j) 

ds  ds  ds 

dW  dO  Trd0  TT  .  nd(j) 

- = - V  —  + 17 sin#  — 

ds  dt  ds  ds 

Just  as  in  the  first  three  equations  of  motion,  these  equations  will  be  useful  later  in  the 
method  of  characteristics.  All  partial  derivatives  here  are  also  with  respect  to  either  time 
or  arc  length  (5),  where  5  is  positive  down  the  towline  towards  the  body. 


1.32 


1.33 

1.34 

1.35 


1.36 

1.37 

1.38 
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3.2.  Applying  the  Method  of  Characteristics 

Mathematical  Character  of  Equations. 

A  few  comments  should  be  made  here  regarding  the  behavior  of  parabolic  versus 
hyperbolic  equations  using  a  stepping  method  such  as  the  method  of  characteristics.  Of 
most  importance,  perhaps,  is  the  fact  that  only  hyperbolic  PDE’s  can  be  used  with  the 
method  of  characteristics.  Tannehill,  Anderson,  and  Pletcher  (Tannehill,  Anderson,  and 
Pletcher,  1997:26-32)  discuss  the  behaviors  of  hyperbolic  versus  parabolic  PDE’s,  noting 
the  dependence  that  each  equation  has  with  respect  to  a  length-time  mesh.  They  note  that 
hyperbolic  equations  are  dependent  only  upon  the  solution  contained  within  the  bounds 
set  by  the  characteristics  at  the  previous  timestep.  This  solution  is  represented  in  Figure 
3.2-1  (similar  figure  in  Tannehill,  Anderson,  and  Pletcher,  1997:27),  where  it  can  be  seen 
that  a  single  solution  exists  at  the  intersection  of  the  characteristic  lines. 
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Figure  3.2-1  Hyperbolic  Characteristic  Lines  &  Domain  of  Interest 
Tannehill,  Anderson,  and  Pletcher  note,  however,  that  parabolic  equations,  unlike 
hyperbolic  equations,  are  dependent  upon  the  entire  physical  domain,  and  thus  have  no 
characteristic  lines;  or  rather,  the  slopes  of  the  characteristic  lines  are  zero.  This  is  why 
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the  method  of  characteristics  cannot  be  applied  to  parabolic  PDE’s.  The  dependence  of  a 
parabolic  solution  is  shown  in  Figure  3.2-2  (similar  figure  in  Tannehill,  Anderson,  and 
Pletcher,  1997:3 1),  where  it  can  be  seen  that  an  infinite  number  of  solutions  exist  at  ti\  or 
rather,  the  solution  is  dependent  on  everything  at  its  current  timestep. 
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Figure  3.2-2  Parabolic  Domain  of  Interest 


In  order  to  detennine  whether  an  equation  is  hyperbolic  or  parabolic,  an  attempt 
must  be  made  to  find  the  characteristic  roots  (these  are  the  slopes  of  the  characteristic 
lines).  If  the  roots  exist  and  are  distinct  and  real,  the  equation  is  hyperbolic.  If  there  is 
only  one  distinct  real  root,  the  equation  is  parabolic.  Another  case  exists  where  the  roots 
are  complex  or  no  real  root  exists,  causing  the  equation  to  by  elliptical.  This  last  case 
turns  out  to  not  be  relevant  here,  however. 

Equations  1.28  and  1.37  are  both  parabolic  (Schram  and  Reyle,  1968:  216), 
whereas  the  rest  of  the  equations  of  motion  are  hyperbolic.  The  zero  characteristic  roots 
that  appear  in  the  parabolic  equations  are  due  to  the  towline’s  inextensible  assumption. 
Equations  1.28  and  1.37  can  be  shown  to  be  parabolic,  with  the  rest  being  hyperbolic,  by 
a  little  work  and  inspection. 
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Tannehill,  Anderson,  and  Pletcher  note  the  forms  of  hyperbolic,  parabolic,  and 
elliptical  PDE’s  (Tannehill,  Anderson,  and  Pletcher,  1997:25)  in  the  following  manner. 
Hyperbolic  PDE’s  exist  in  the  two  fonns  of  Equations  2.1  and  2.2: 


dlLJlL=Q 

dt 2  ds2  Sl 


¥  ¥ 
dt’  ds’ 


2.1 


where  g(  )  denotes  some  function  of  the  included  variables. 

To  determine  if  Equations  1.27  through  1.29,  and  Equations  1.36  through  1.38  are 
hyperbolic,  parabolic,  or  elliptical,  they  need  to  be  converted  into  a  second  order  fonn. 

This  is  done  by  noting  that  the  velocity  terms  of  the  equations  in  question  can  be  replaced 

2  2 

with  their  respective  position  derivative  (i.e.,  El  =  dX/  dt,  thus  dU/dt  =  d  XI  d  t).  The 
angles  can  also  be  converted  into  position  terms  by  manipulating  Equations  1.12  through 
1.14  into  the  following: 

^tan-'lA)  2.6 
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1 


dz 


2.7 


From  here,  one  can  convert  all  six  equations  of  motion  to  fonns  that  have 
derivatives  of  solely  position  variables.  Another  transformation  between  space  and 
towline  fonns  would  be  required  for  consistency,  although  to  simply  determine  the 
character  of  the  equations  by  inspection  this  would  not  be  necessary  since  it  would  not 
add  more  derivative  terms.  In  other  words,  x  and  X are  interchangeable  for  the  inspection 
technique. 

The  details  of  this  procedure  are  not  vital  to  understanding  the  methodology  of 
this  paper,  and  will  thus  not  be  included  here.  Suffice  it  to  say  that  Equations  1.28  and 
1.37  are  both  parabolic,  whereas  the  rest  of  the  equations  of  motion  are  hyperbolic.  Since 
only  hyperbolic  equations  can  be  applied  to  the  method  of  characteristics,  the  parabolic 
equations  cannot  be  transformed  in  the  subsequent  sections. 

Finding  Roots  of  Equations. 

The  roots  of  the  hyperbolic  equations  can  be  found  by  setting  up  a  generalized 
Eigenvalue  problem.  Generalized  Eigenvalues  are  defined  as 

[u][M]  =  A[u][JV]  2.8 

where 

det(M  -  XN)  =  0  2.9 

[ M]  and  [A]  are  coefficient  matrices  of  a  system  and  [t>]  is  chosen  such  that  equation  2.8 
holds  true  (in  reality,  [  u  ]  is  said  to  be  non-zero  for  the  system  to  have  Eigenvalues,  thus 
its  values  do  not  matter  for  finding  roots,  since  it  follows  that  the  determinant  must  be 
zero  for  the  system  to  hold).  Solving  for  X  finds  the  roots  of  the  system.  Schram  and 
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Reyle  model  the  system  of  hyperbolic  equations  based  on  a  technique  they  developed 
from  Ames  (Schram  and  Reyle,  1968:216-217).  In  this  they  converted  the  four 
hyperbolic  equations  (1.27,  1.29,  1.36,  1.38)  into  one  single  system  of  matrices. 

Mfl+fiM+RHO 


where 


[M] 


0  0 

0  0 

1  0 
0  1 


Tcos9 

M 

0  -- 

A 

-(V  cos  9  -  W  sin  9)  0 

-U  sin  9  V 


-1  0  (V  cos  9  -  W  sin  9)  0 

0  -1  U  sin#  -V 

0  0  —  cos  9  0 

0  0  0  1 


M 


u 

w 

<t> 

9 


IP] 


1 


Fx  +  Wt  cos  (f) 
Wt  sin  (j)  sin  9  +  Fz 
0 
0 


2.11 


2.12 


2.13 


2.14 


and  the  5  and  t  subscripts  refer  to  the  partial  derivatives  of  the  matrix  variables. 

From  this,  one  can  find  the  general  Eigenvalues  of  the  system,  which  represent 
the  characteristic  roots.  This  is  done  by  solving  for  X  in  equation  2.9.  The  roots  of  the 
system  are  thus  represented  by 
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which  are  each  repeated  twice  (four  equations  gives  four  roots).  The  subscripts  a  and  P 
refer  to  the  characteristic  directions.  These  values  turn  out  to  be  the  speed  at  which 
disturbances  are  propagated  along  the  towline.  In  other  words, 


This  is  consistent  with  the  value  that  can  be  found  in  any  textbook  for  propagation  of 
disturbances  along  an  inextensible  line  in  tension. 


Finding  Characteristic  Equations. 

Schram  and  Reyle  (Schram  and  Reyle,  1968:217)  note  that  the  characteristic 
equations  can  be  determined  by  some  manipulation  of  the  hyperbolic  equations.  This  can 
be  done  by  multiplying  Equations  1.27  and  1.29  by  dt  and  adding  them  to  their 
counterpart  Equations  1.36  and  1.38  multiplied  by  ds.  This  provides  two  equations  that 
have  partial  derivative  terms  with  respect  to  both  s  and  t.  By  noting  that  the  total 
derivative  is  defined  as 

df  (s,t)  =  —ds  +  —  dt  2.18 

ds  dt 

the  partial  derivative  terms  of  these  two  equations  can  be  combined  (this  process  involves 
some  simple  but  lengthy  equations  and  will  thus  not  be  included  here).  Replacing  the 
ds/dt  term  with  the  characteristic  root  of  Equation  2.17,  this  provides  two  equations  with 
no  partial  derivatives: 
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dU  +  d<f> 


cos  6  -  V  cos  6  -  W  sin  6 


— (Fx  +  Wt  cos  </>)  =  0 
M 


2.19 


dW  +  dd 


U sin Od(f>-  —  {Wl  sin^sin#  +  Fz)  =  0 


2.20 


These  equations  are  the  total  derivatives  in  the  characteristic  directions.  It  should  be 
noted  that  the  parabolic  equations,  had  they  been  attempted  to  be  converted  in  this 
manner,  would  have  remaining  partial  derivative  terms  and  would  thus  have  no  benefit  in 
converting  their  form. 

By  taking  the  derivative  of  these  equations  with  respect  to  the  characteristic 
directions,  one  finds  the  characteristic  equations  that  can  be  used  for  integration: 

<^_  +  ^0_\Fa  C0SQ-y  cos#-IFsin#]--^-  — (Fv  +  Wt  cos (j))  =  0 
da  da  da  // 


—  +  ^\  Fr  cosd-V  cos  6-  W  sin  0~\ 
dp  dp1  p  J 


dW  d0  r  ,  d(j)  dt  1 

- 1 - \V  -Fa\-U  sinfe1^ - 

da  da  da  da  n 


—(Fx+Wt  cos  (f>)  =  0 
dp  // 

(Wt  sin^sin @  +  Fz)  =  0 


—  +  —  \V-Fp-\-Usm9^-—-(Wtsm(j)smO  +  Fz)  =  0 
dp  dp1  dp  dp  ii  ‘ 


It  should  be  noted  that  a  and  P  are  the  characteristic  directions  in  the  z-y  plane,  but  since 


2.21 

2.22 

2.23 

2.24 


the  characteristic  roots  are  repeated,  these  changes  are  the  same  in  both  the  z-y  and  the  x- 


y  planes  (Schram  and  Reyle,  1968:218). 
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3.3.  Numerical  Methods 


Setting  up  the  Mesh  and  the  Characteristic  Equations . 

As  mentioned  in  Section  2.6,  the  numerical  method  used  in  this  paper  is  one 
developed  by  Courant,  Isaacson,  and  Rees,  presented  in  Ames  (Ames,  1965:  446),  and 
applied  to  this  type  of  problem  by  Schram  and  Reyle  (Schram  and  Reyle,  1968:218). 

This  technique  uses  a  preset  space-time  mesh  in  order  to  reduce  the  amount  of 
calculations  required.  As  a  result,  the  mesh  size  must  be  sufficiently  small  enough  to 
ensure  straight  characteristic  lines  between  timesteps. 

Using  this  understanding  as  background,  one  may  rewrite  Equations  2.21  through 
2.24  into  a  set  of  simplified  characteristic  ODE’s 


dU  +  Gad(j)  +  Hadt  =  0 
dU  +  Gpd(j)  +  H pdt  —  0 
dW  +  Jad6  +  Kad(j)  +  Ladt  =  0 
dW  +  J pdQ  +  Kpd(f>  +  Lpdt  —  0 

where  the  coefficients  are  defined  as 

Ga  =W sind-V cos 6  +  Fa  cos 6 

Gp  =  W  sin  6  -  V  cos  0  +  Fp  cos  0 


J „  =  V  —  F„ 


J  =V-F 

J  p  v  1  p 


Ka  =  —U  sin  0 


Kp  =  -U  sin  6 


3.1 

3.2 

3.3 

3.4 


3.5 

3.6 

3.7 

3.8 

3.9 
3.10 
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3.11 


Ha=-(~FX-Wt  cos^) 
M 


La  =  —  (  -Fz  -  Wt  sin  0  sin  6?) 


For  the  towline,  the  following  coefficients  are  defined  as 

Hp  =  —{-Fx  -Wt  cos<ft) 

A 


Lp  =  —  (~FZ  -W,sm0  sin  9) 
A l 


At  the  body,  however,  these  are  defined  as 


HB=—(-Fx-Wbcosd>) 
p  Mb 


Lp  =  -ffbsin^sin#) 


where  the  forces  are  calculated  on  the  body,  noting  that  only  the  /3  characteristic  lines 
affect  the  lower  boundary,  and  of  these  /3  characteristic  parameters,  only  the  force  terms 
become  different  (thus  only  two  terms  need  to  be  changed).  It  should  also  be  noted  that 
the  characteristic  root  changes  at  the  lower  boundary,  thus, 


Fa 


3.12 


3.13 

3.14 


3.15 

3.16 


3.17 


at  the  towed  body. 

These  equations  are  applied  to  every  point  down  the  towline,  in  order  to  “march” 
the  data  to  the  next  timestep.  This  is  represented  in  Figure  3.3-1.  Points  a,  b,  and  c  are 
points  along  the  towline  where  data  is  known.  Points  p  and  q  are  set  by  the  Fa  and  Fp 
characteristic  lines  by  noting  that  these  lines  are  defined  as  dsldt.  Thus,  ds,  which  is  the 
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distance  along  the  5-axis  between  both  q  and  b  as  well  as  b  and  p,  can  be  defined  as  either 
\Fa\dt  or  \Ffj\dt. 


Figure  3.3-1  Character  Mesh  with  Characteristic  Lines 
The  data  associated  with  the  points  p  and  q  can  be  found  by  interpolating  between 
the  known  data  points  (a,  b,  and  c)  in  the  following  manner: 

f(q)  =  -1- [f(a)ds  +  f(b)(As  -  ds)] 

As 


3.18 


f(p)  =  ^7 [f(c)ds  +  f(b)(As  -  <&)] 


3.19 


where^  )  refers  to  any  type  of  data  stored  at  the  respective  points  (i.e.,  velocity). 

As  can  be  seen  in  Figure  3.3-1,  the  stability  of  the  preset  mesh  is  restricted  by  the 
points  p  and  q  being  between  a  and  c.  The  maximum  values  of  the  characteristic  lines  are 
thus  limited  in  this  mesh  to  At/ As.  These  values  force  q  to  line  up  with  a,  and p  to  line  up 
with  c.  In  other  words,  the  stability  of  this  system  is  restricted  by 


3.20 


3.21 


42 


which  forces  ds  to  not  go  beyond  As.  Thus  the  mesh  remains  stable. 

This  setup  produces  two  sets  of  equations  (one  from  each  characteristic  line)  for 
every  point  in  time  and  space  for  all  data  (i.e.,  velocity)  except  at  the  boundary 
conditions,  which  only  have  one  characteristic  line  at  each  end  (the  upper  boundary 
condition  only  uses  Fa  and  the  lower  boundary  condition  only  uses  Fp).  By  going  back 
and  writing  the  characteristic  equations  (Equations  3.1  through  3.4)  in  difference  notation 
form,  the  equations  for  numerical  integration  along  the  line  are  completed.  It  should  be 
noted  that  all  the  coefficient  variables  (Equations  3.5  through  3.16)  are  evaluated  at  point 
b  on  the  mesh,  which  Ames  notes  as  well  (Ames,  1965:447). 

Characteristic  Equations: 


Ur-Up+GM-tP)  +  H„dt  =  0 

3.22 

Ur-Uq+GfiWr-tq)  +  Hfidt  =  0 

3.23 

Wr-Wp+Ja{Or-Op)  +  KaUr-tp)  +  Ladt  =  0 

3.24 

K  -  wq  +  Jp  ( er  -  eq ) + kp  ( cj>r  - 1 ) + lp  dt  =  o 

3.25 

These  equations  can  be  combined  in  order  to  solve  for  the  four  unknowns:  Ur,  Wr,  </> 
and  0r. 

Since  the  two  parabolic  equations  (Equations  1.28  and  1.37)  were  not  transformed 
by  the  method  of  characteristics,  these  must  be  directly  translated  here  into  difference 
notation  so  that  they  can  be  used.  Schram  and  Reyle  use  a  central  differencing  technique 
in  order  to  ensure  a  more  accurate  result  (Schram  and  Reyle,  1968:219),  nothing  that 
Equation  3.26  is  a  central  differencing  about  an  averaged  value  (halfway  between  e  and 
r). 
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Parabolic  Equations: 


,  W  +  W 
V  -  V  =  — _ — 

2  j 


e  r 


(9.-9,)- 


fU,+U,'' 

V  2  , 


COS 


C 0+6,  > 


V  £  ) 


(k-h) 


Td~Te=Ss 


£_ 

St 


[{K-Vb)~Wr  iQr  ~eb)  +  Ur  C0S  6r  (</>r  ~  <!>b  )]  ~  FY,r  +  W,  sil1  <Pr  C0S  dr 


where  Ss  is  the  distance  between  the  two  points,  which  becomes  2A.v  here,  and  St  is  the 
distance  between  time  steps,  which  becomes  At. 

These  equations  must  be  applied  after  the  characteristic  equations  have  been  used 
to  find  the  values  down  the  entire  line.  The  boundary  conditions  (discussed  later)  must 
be  defined  at  the  upper  boundary  in  advance  for  Equation  3.26.  This  equation  can  then 
be  used  for  all  values  including  the  lower  boundary  (a  different  equation  will  be  used 
later,  however).  The  boundary  conditions  must  also  be  defined  at  the  lower  end  in 
advance  for  T  since  it  is  to  be  found  starting  at  the  towed  body  and  calculated  back  up  the 
towline  to  the  aircraft  (see  Figure  3.5-1  at  the  end  of  this  chapter).  Since  the  central 
differencing  technique  of  Equation  3.28  cannot  be  applied  at  the  first  point  up  the  towline 
from  the  lower  boundary,  a  forward  differencing  scheme  is  used  instead  (Schram  and 
Reyle,  1968:219). 


Td~TL=  Ss 


Yt  [iVd  -  va ) -  wd (°d  - 6 a )  +  ^ Ud  cos ed {(f>d-(f>a)\- FYd  +  Wt  sin (f)d  cos Qd 


where  Ss  is  A.v  and  At  is  At. 

This  produces  a  set  of  six  equations  (plus  a  seventh,  but  only  due  to  the  issue  of 
central  differencing)  that  can  describe  the  behavior  down  the  entire  line  except  for  at  the 
boundaries.  The  boundary  conditions  for  the  mesh  used  here  exist  at  5  =  .sa/c  and  5  =  Sbody. 


3.26 

3.27 


3.28 
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In  addition,  the  initial  conditions  at  t  =  0  must  be  known  before  any  kind  of  stepping 
technique  is  used. 

Initial  Conditions. 

The  initial  conditions  are  detennined  by  prior  research  using  MATLAB®  code  to 
develop  a  steady  state  solution  (Richardson,  2005:Appendix  A).  The  method  by  which 
this  code  was  developed  is  discussed  in  Section  2.1.  Another  way  to  set  the  initial 
conditions  is  by  setting  the  angles  of  the  towline  and  its  velocities  in  all  three  axes 
manually.  The  tension  can  be  found  up  the  entire  line  by  calculating  the  body  and 
towline  forces  (shown  later)  from  the  known  angles  and  velocities. 

Boundary  Conditions. 

In  order  to  find  the  values  at  the  upper  and  lower  boundaries,  new  equations  must 
be  developed  to  work  with  the  characteristic  lines.  This  is  due  to  the  fact  that  there  must 
be  two  equations  for  every  point  in  order  to  detennine  the  values  at  that  point.  The 
characteristic  lines  carry  data  to  the  boundary  conditions  as  can  be  seen  for  the  upper 
boundary  conditions  in  Figure  3.3-2. 
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Figure  3.3-2  Characteristic  Line  to  find  upper  BC  data 

Since  it  is  assumed  that  the  towline  at  the  point  where  it  connects  to  the  aircraft 
has  the  same  motion  as  the  aircraft,  the  new  equations  for  the  upper  boundary  conditions 
in  towline  form  become  the  transformation  of  coordinates  from  Equation  1.6.  Thus,  the 
upper  boundary  conditions  can  be  described  by: 


V 


=  [A]  M„ 


3.29 


where 


refers  to  the  velocity  matrix  in  towline  coordinates  at  the  upper  limit  (i.e., 


the  aircraft),  and  [v]a  refers  to  the  velocity  matrix  in  spatial  coordinates  at  the  upper 

limit.  By  multiplying  through,  the  upper  boundary  equations  are  found: 

Upper  Boundary  Conditions: 

Uu  =uusm0"-vucos0u  3.30 

Vu  =  uu  cos  6U  cos  </>u  +  vu  cos  0u  sin  <f>u  +  wu  sin  9U  3.31 

Wu  =  -uu  sin  0U  cos  (f)u  -  vu  sin  6U  sin  (f)u  +  wu  cos  6u  3.32 

Equation  3.22 
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Equation  3.24 


The  equations  are  all  combined  to  find  the  values  at  the  upper  boundary:  Uu,  Vu, 
W u,  <f>  u,  and  0U  (Tu  will  be  found  later  by  3.27).  Schram  and  Reyle  (Schram  and  Reyle, 
1968:219)  note  that  an  iterative  procedure  must  be  used  to  find  (j>  u  and  0„  such  that  these 
equations  are  satisfied. 

The  lower  boundary  conditions  become  a  little  trickier.  However,  a  few 
assumptions  can  be  made  that  simplify  the  process.  The  characteristic  line  used  at  the 
lower  boundary  can  be  seen  in  Figure  3.3-3. 


t 


Figure  3.3-3  Characteristic  Line  to  find  lower  BC  data 
With  the  assumption  that  the  towline  at  the  towed  body  has  the  same  motion  as 

the  body,  one  can  develop  the  lower  boundary  conditions.  This  requires  a  little 
understanding  of  the  geometry.  For  the  sake  of  simplicity,  the  body  will  be  assumed  to 
be  rigid,  small,  and  always  parallel  to  the  freestream  direction,  thus  producing  negligible 
moment  effects.  In  other  words,  the  body  will  be  assumed  to  be  reducible  to  a  point  mass 
in  how  its  behavior  transfers  to  the  towline. 
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By  taking  Equations  1.36,  1.37,  and  1.38  and  writing  them  in  finite  difference 
form  using  a  central  differencing  technique  about  point  d  in  Figure  3.3-3,  one  finds  three 
equations  of  motion  at  the  lower  boundary: 

Ss 

UL  =U  +—(</>L  -</>b)  cos  ed  +  (Vd  cos  ed  -  Wd  sin  6d  ){(/>L  -  <j>  ) 

St 


vL  =  vg-  ud  cos  ed  (A -0g)+wd  (0L  -  eg ) 


S  S 

trL=Wg-  —  (0L-0b)-Vd(0L-0)  +  Udsin0M-0) 


5t 


3.33 

3.34 

3.35 


Central  differencing  is  used  due  to  the  fact  that  forward  differencing  creates  a  set  of  five 
nonlinear  equations  that  have  no  exact  symbolic  solution.  Central  differencing  prevents 
the  need  for  using  iteration  to  solve  the  system  for  a  best  value.  It  should  be  noted  that 

Ss  As 

due  to  this  setup,  Ss  is  set  to  twice  the  As  value.  Thus,  —  =  2 —  . 

St  At 

Combining  these  three  equations  with  Equations  3.23  and  3.25,  five  equations 
have  now  been  detennined  for  the  lower  boundary,  which  contain  five  unknown 
variables,  (f>,  can  be  determined  explicitly  and  9l  becomes  a  function  of  only  <j)L  from 

these  equations  by  first  creating  two  new  equations.  This  is  found  by  rewriting  Equations 
3.23  and  3.25  in  difference  notation  at  the  lower  boundary: 

UL=Uq-GM-tP)~Hadt  3.36 

WL  =  wq  -  Jp  ( 9 L  -  9 p  ) -  Kp  (4  -  <t>q )  ■ -  L,dt  3.37 


These  two  equations  are  then  combined  with  Equations  3.33  and  3.35  respectively  to 
form  two  new  equations,  the  first  dependent  solely  on  (j>L  and  the  second  on  both  9l  and 


<t>L- 
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3.38 


Ss 

ug  + — (fa  -  fa ) cos  fa,  +  (K,  cos  ed  -  wd  sin ed  -  </>  ) 


St 


-Ua+GflWL-ta)  +  Hedt  =  0 


Ss 

wg-—{eL-eb)-vd{eL-eg)+ud^eML-(i>s) 

~Wq+  Jp{0L  -6q)  +  Kp(<j)L  -<pq)  +  Lpdt  =  0 


3.39 


As  can  be  seen,  the  only  unknown  in  Equation  3.38  is  (f>L .  Thus,  the  equation  can 


be  rearranged  to  reveal 
Ss 


fa  = 


St 


fa  cos  °j  +(vd  cos  °j  ~  wc  sin  °,i  )fa -Ug-HBdt  +  Ua+ 


9  1 


Ss 

Gp  +  vd  cos  ed  -  Wd  sin  0d  +  —  cos  6d 


3.40 


Similarly,  Equation  3.39  can  be  rearranged  to  reveal 


0L  = 


Ss 

Wg  +  Jtfa  +  VA+Ud  si  nOd(fa  -fa)-Wq  -Jpeq  +  Kp(<fiL -fa)  +  Lpdt 

—  +  y  -j 
St  d  p 


3.41 


By  solving  for  these  two  values,  the  values  for  UL,  VL,  and  WL  can  now  be  solved. 

Thus,  the  lower  boundary  conditions  are  defined  by  the  following  equations, 
noting  that  values  for  Up  and  W/  can  be  solved  by  two  different  equations  and  that 


Equations  3.40  and  3.41  must  be  solved  first: 

Lower  Boundary  Conditions: 

Equation  3.40 
Equation  3.41 
Equation  3.36  or  3.33 
Equation  3.37  or  3.35 
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Equation  3.34 


The  values  for  Ul,  Vl,  W l,  <t>  l,  and  dL  have  been  detennined.  These  new  values  are, 
then,  used  to  solve  for  tension. 

Body  Forces. 

The  forces  on  the  body  must  be  found  in  order  to  model  the  towline’s  movement 
and  tension.  This  is  best  done  by  finding  a  coordinate  system  aligned  with  the  body. 
Since  the  body  is  assumed  to  be  attached  to  the  towline  such  that  it  has  no  rotation,  the 
body  axes  become  the  towline  axes  at  attachment.  Schram  developed  a  manner  in  which 
to  analyze  rotation  in  the  body  (Schram,  1968:1 19-120)  using  Euler  angles,  which  will 
not  be  used  or  discussed  here.  His  method  becomes  quite  rigorous,  and  the  assumptions 
here  involve  negligible  moment  effects  (thus  no  angular  acceleration  is  included).  A 
manner  in  which  to  account  for  angular  acceleration  is  included  in  Appendix  A.  The 
body  orientation  is  represented  in  Figure  3.3-4  and  Figure  3.3-5.  Again,  it  should  be 
noted  that  6  is  in  a  unique  plane  set  by  (j) ,  and  could  be  coming  in  or  out  of  the  page. 
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To  calculate  the  forces  in  each  axis,  the  drag  components  must  first  be  found. 
Since  it  is  assumed  that  the  body  is  always  aligned  with  the  freestream  flow,  drag  is 
always  acting  along  the  body,  which  allows  for  the  use  of  a  single  drag  coefficient.  The 
values  will  be  calculated  in  the  spatial  coordinate  system,  since  weight  will  have  to  be 
added  later.  The  spatial  velocities  can  be  found  from  the  towline  velocities  at  the  lower 
boundary  by  the  transformation  of  Equation  1.7.  Noting  that  the  total  drag  force  as  a 
function  of  components  is  defined  as  DTOt  =  Dx  +  Dy  +  Dz2 ,  and  velocity  must  be 
decomposed  based  on  these  forces, 


Dx  ^  P^LUS'D,B 

Dy  =  jP^LVLCD,B 

Dz  =  ~P^LWL^D,B 


-F 


-F, 


-F_ 


3.42 

3.43 

3.44 


where 


Vl=Fl1+vl1+W  3.45 

Note  that  the  negative  sign  on  the  force  tenns  is  due  to  their  orientation  being  in  the 
positive  direction  in  each  axis.  Since  drag  acts  the  same  direction  as  airflow,  these  force 
tenns  must  have  a  negative  sign. 

It  is  assumed  that  there  is  no  lift  on  the  body.  A  lift  force  can  be  easily  added  into 
the  equations,  although  it  will  not  be  done  here.  The  aerodynamic  force  tenns  are 
therefore  the  drag  terms  and  were  put  into  the  above  equations  for  consistency  with  prior 
equations  (3.15,  3.16). 
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Additional  forces  based  on  the  body’s  acceleration  must  by  found  by  Newton’s 


Second  Law  (Equation  1.8).  The  acceleration  can  be  found  through  the  use  of  a 
backward  differencing  method  outlined  in  Tannehill,  Anderson,  and  Pletcher  (Tannehill, 
Anderson,  and  Pletcher,  1997:50) 

3M(»-4M(>-l)  +  «(>-2)  +  q  2 
2  dt 


as  well  as  a  simple  backward  differencing  method 

u(t)-u(t-  1) 

u  = - 

dt 


+  0{dt ) 


3.47 


for  the  first  timestep.  The  more  complex  method  is  used  in  order  to  reduce  truncation 
errors,  which  are  one  order  higher  than  the  simpler  method. 

Thus,  the  force  tenns  due  to  inertia  for  all  three  axes  become: 

Faccel,x  =  ilIMh 
Faccel,y  =  T LM^ 

Faccel,z  = 


3.48 

3.49 

3.50 


It  should  be  noted  that,  due  to  Newton’s  Third  Law,  these  forces  act  away  from  the 
aircraft.  In  other  words,  a  positive  acceleration  creates  a  negative  force.  Also  of 
important  mention  is  the  fact  there  can  be  no  compressive  forces  in  the  line.  Thus,  any 
acceleration  that  would  theoretically  compress  the  line  cannot  be  factored  into  the  forces 
for  tension.  Since  compression  is  impossible,  the  acceleration  has  no  effect  on  tension. 
This  produces  abnormal  results  in  the  code,  thus  the  mass  and  drag  coefficient  must  be 
monitored. 

These  forces  can  all  be  summed  up  to  detennine  the  total  forces  on  the  body  (not 
including  tension),  noting  that  the  drag  acts  in  the  same  direction  as  airflow,  the 
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acceleration  is  oriented  away  from  the  aircraft,  and  the  force  due  to  weight  must  be  added 
to  the  v-axis  tenn: 

TP,  =  -A  LMb  -  i  pfLuLCD,  ^  3.51 

TFr=-vLMb-^pfLvLCD,i^-Wb  3.52 

TF,  =  -wLMb-^pVLwLCDti^-  3.53 

The  tension  acts  purely  in  the  towline  T-axis,  and  since  the  angles  were  found  for 
the  lower  boundary,  the  transformation  of  Equation  1.6  can  be  applied.  This  yields: 

TL  =  -{TFx  cos </> cos 6 +  TFy  sin^cos^  +  777.  sin#)  3.54 

where  the  negative  sign  is  due  to  the  tension  being  opposite  the  direction  of  total  force 
such  that  the  sum  of  all  forces  in  the  T-axis  is  zero.  Thus,  TL  is  found  at  the  lower  end, 
and  can  be  calculated  back  up  the  towline  by  the  use  of  Equations  3.27  and  3.28,  noting 
the  towline  forces  derived  in  the  next  section. 

Towline  Forces. 

The  towline  forces  (circular  towline  with  no  rotation)  can  be  calculated  in  a 
similar  manner.  The  drag  forces  on  the  towline  are  modeled  in  Figure  3.3-7  and  Figure 
3.5-1. 
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Figure  3.3-6  Towline  Force  Components  (x-y  plane) 


y 


Figure  3.3-7  Towline  Force  Components  (z-y  plane) 

One  major  difference,  however,  is  that  the  towline  is  not  considered  to  be  facing 
the  freestream  velocity.  Although  drag  is  to  be  considered  a  total  force  (i.e.,  it  must  first 
be  detennined  based  on  total  velocity  before  decomposition  into  axial  components), 
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Note  again,  the  force  terms  are  to  be  in  the  positive  direction,  hence  their  negative  sign. 
Since  there  are  no  lift  forces,  the  aerodynamic  forces  become  the  drag  forces,  and  can  be 
applied  to  the  previous  equations  (3.11  through  3.14,  3.27,  and  3.28). 
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3.4.  Heat  Transfer 


The  entire  system  is  assumed  to  be  in  a  predetennined  temperature  field.  Thus, 
the  temperature  of  the  air  is  known  at  every  point  in  space  and  time.  In  order  to 
determine  the  heat  transfer  in  the  system,  an  energy  equation  must  be  found  that  relates 
the  rate  of  internal  energy  transfer  between  the  air  and  the  towline.  The  orientation  of  the 
control  volume  at  each  incremental  length  of  towline  is  represented  in  Figure  3.4-1. 


ds 


Figure  3.4-1  Airflow  Over  Towline 

Due  to  conservation  of  energy,  the  energy  balance  equation  may  be  written  as 

E  - E  +  E  —  E  =0 

in  out  generated  stored 

where  all  values  are  in  terms  of  the  rate  of  internal  energy.  Since  the  towline  is  not 
generating  any  energy,  this  equation  becomes 

E  =  E 

in  stored 

assuming  there  is  one  dimensional  heat  transfer  from  the  air  to  the  towline.  If  it  were  the 
other  way  around,  the  sign  on  Ein  would  be  switched  (essentially  becoming  -Eout ),  so 

Equation  4.2  holds  true  for  both  heat  being  transferred  to  the  line  and  away  from  it. 

To  solve  for  the  energy  rate  into  the  towline,  Newton’s  law  of  cooling  can  be 

used: 
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E,n=hAs(  T.-T,) 


4.3 


where  As  is  the  total  surface  area  over  which  the  air  flows,  and  h  is  the  convection 
coefficient.  is  used  here  since  the  heat  transfer  is  between  the  far  field  flow  and  the 

towline  surface  ( Ts ). 

Since  a  thin  boundary  layer  forms  about  the  towline,  film  temperature  will  be 
used  to  find  the  properties  of  the  air  at  the  towline’s  surface.  Film  temperature  is 
assumed  to  be  the  average  temperature  between  the  far  field  and  the  towline  and  is 
defined  as 


T  -T 

rp  _  00 _ S_ 


4.4 


The  convection  coefficient  can  be  defined  as 


h=^-k 


D 


TL 


where  Nit  is  the  Nusselt  number,  k  is  the  thermal  conductivity  of  the  air,  and  Dtl  is  the 
total  distance  across  the  towline  (see  Figure  3.4-1).  The  distance  across  the  towline  can 
be  defined  in  the  following  manner,  noting  the  perpendicular  and  parallel  components  of 
the  velocity  across  the  line,  which  create  an  angle  at  which  the  airflow  acts: 


V. 


perpendicular 


=  ^U2+W2 


4.5 


4.6 


v  =  v 

parallel 


4.7 


FlowAngle  =  tan 


v 

perpendicular 

.  V 

\  parallel  J 


4.8 


Thus, 
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4.9 


Dtl  = 


sin  {FlowAngle) 


for  airflow  that  completely  crosses  the  towline  in  the  perpendicular  direction  (i.e.,  it  does 
not  completely  cross  the  parallel  direction).  Otherwise, 

ds 


Dtl  = 


cos  {FlowAngle) 


for  airflow  that  completely  crosses  the  towline  in  the  parallel  direction.  The  towline’s 
true  diameter  ( dTL  )  will  be  used  for  now,  however,  as  will  be  shown  later  in  the 
discussion  of  the  Nusselt  number. 

The  value  for  k  can  be  calculated  directly  from  the  film  temperature.  This  is 
shown  by  a  power  fit  to  experimental  data  (personal  communication,  Ralph  Anthenien) 


as 


k  =  (0.0002235)T 


0.8302 


4.10 


4.11 


where  T,  is  in  Kelvin. 

The  Nusselt  number  is  calculated  experimentally.  Currently,  an  approximation 
for  the  Nusselt  number  could  not  be  found  for  a  cylinder  with  different  angles  to  the 
freestream.  Also,  the  boundary  layer  cannot  be  assumed  to  be  very  small  down  the 
towline  (boundary  layer  thickness  increases  along  the  towline),  which  will  affect  heat 
transfer.  Thus,  only  the  perpendicular  component  of  heat  transfer  will  be  calculated  as  if 
the  only  flow  across  the  towline  is  perpendicular  (i.e.,  thin  boundary  layer  assumption), 
leaving  a  better  Nusselt  number  approximation  for  future  work.  This  is  consistent  with 
the  rest  of  the  paper,  which  assumes  no  drag  parallel  to  the  towline.  Thus,  the  Dtl  term 
in  Equation  4.5  becomes  dn-  For  perpendicular  flow  across  a  cylinder,  Nu  for  Reynolds 
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numbers  within  the  range  1  <  ReD  <  104  is  shown  by  Kramer  as  reproduced  in  White 
(White,  1991:186)  as 

Nu  =  0.42/V02  +  0.57  Prm^R^ 

This  fonnula  was  detennined  from  a  curve  fit  of  experimental  data.  Pr  is  the 
Prandtl  number,  assumed  to  be  ~0.7  for  air  in  the  range  of  -275-850  Kelvin.  The  Prandtl 
number  changes  very  little  as  a  function  of  temperature,  and  this  approximation  could  be 
applied  to  even  a  greater  range  of  temperatures.  The  Reynolds  number  is  calculated  by 


ReD 


Vperpendiculard 


TL 


V 

where  v  is  the  kinematic  viscosity,  which  can  calculated  directly  by  a  fit  to  experimental 
data  (personal  communication,  Ralph  Anthenien)  based  on  the  film  temperature: 


t/  =  8x  10-10T/'7235 

where  Tf  is  again  in  Kelvin.  Future  development  of  a  Nusselt  number  approximation 

that  accounts  for  parallel  flow  will  also  have  to  develop  a  new  Reynolds  number 
approximation  as  well,  assuming  that  Nu=J[Re). 

In  order  to  detennine  the  stored  energy  in  the  towline,  specific  heat  is  used. 
Specific  heat  is  defined  as  the  amount  of  heat  per  unit  mass  required  to  raise  the 
temperature  one  degree  Kelvin.  Since  heat  can  be  defined  by  the  change  in  stored 
internal  energy  over  the  change  in  temperature,  specific  heat  for  the  towline  is 
represented  as 


c  = 


J stored 


pVol 


dl 

dt 


4.12 


4.13 


4.14 


4.15 
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Here,  p  and  V»i  are  the  density  and  volume  of  the  towline  respectively,  and  d  T  /dt  is  the 
change  in  temperature  over  time.  The  dt  tenn  comes  in  because  the  stored  internal 
energy  is  a  rate  as  well. 

Rearranging  this  gives  the  change  in  energy  storage  due  to  the  temperature  change 

4.16 

dt 

where  the  specific  heat  value  used  in  this  analysis  was  an  approximate  value  for  steel  of 
450  J/kg-K.  Combining  Equations  4.2,  4.3,  and  4.16  yields  the  energy  equation  at  the 
surface  of  the  towline 


hAs(Tf-Ts)  =  cpVo,d^ 

Thus,  temperature  can  be  calculated  at  any  point  by  knowing  the  position, 
velocity,  and  previous  temperature  of  the  towline: 


4.17 


T  =  T  +- 

1  s,l+ 1  1  T 


cpVoi 


4.18 


where  Tv  is  the  temperature  at  the  surface  at  the  current  time  and  Tv  Ml  is  the  new 
surface  temperature  at  the  new  time. 

A  new  film  temperature  can  then  be  calculated  based  on  the  new  surface 
temperature  value  through  the  use  of  Equation  4.4.  This  is  used  to  calculate  the  air 
density  in  the  boundary  layer  around  the  towline,  which  is  used  in  turn  to  find  the  forces 
on  the  towline.  Density  can  be  calculated  by  fit  to  experimental  data  in  SI  units  (personal 
communication,  Ralph  Anthenien)  as 

/?  =  (357.88)T/1  0041  4.19 

where  it  should  again  be  noted  that  temperature  is  in  Kelvin. 
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3.5.  Integration  Procedure 


Enough  equations  have  now  been  derived  to  model  the  behavior  of  the  towline 
and  body,  as  well  as  the  heat  transfer  to  the  towline,  at  all  points  in  time  and  space. 

These  equations  were  written  into  MATLAB®  code  following  the  numerical  integration 
procedure  as  can  be  seen  in  Figure  3.5-1.  The  code  is  attached  in  Appendices  A  through 
D. 


1 .  Find  upp  er  b  oundary  c  onditions  (U,V,W,  6  ,<3>)  from 
previous  time  and  current  a/c  velocity  (u,v,w). 

2.  Calculate  <£,11,©,  and  W  down  entire  line  from  previous 
time,  then  find  V  from  these  values. 

3.  Find  ©  explicitly  at  lower  boundary  from  previous  time 
and  values  one  step  up  line,  then  find  ®,U,V,  and  W. 

4.  Calculate  T  at  lower  boundary  from  forces  on  the  body, 
then  find  T  up  entire  line  and  at  upper  boundary. 

5.  Find  x,y,z  spatial  positions  along  line. 

6.  Find  temperature  and  density  at  every  point  along  line. 

7.  Go  to  next  timestep. 


Figure  3.5-1  Numerical  Integration  Procedure 


62 


IV:  Results  and  Discussion 


The  MATLAB®  code  was  written  such  that  a  change  in  initial  conditions  and 
system  values  could  be  easily  altered.  The  behavior  of  the  towline  was  initially  modeled 
with  the  aircraft  under  a  constant  speed  with  the  towline  extended  directly  back  behind 
the  aircraft  at  a  constant  angle  down  the  line  of  0.01  (a  zero  angle  gives  a  negative 
tension  error  almost  immediately).  These  results  were  compared  to  the  steady  state 
values  from  Richardson.  This  same  procedure  was  also  applied  using  the  initial 
conditions  from  Richardson  to  test  for  constant  values  under  steady  state.  The  system 
was  then  modeled  under  different  perturbations,  starting  with  the  constant  angle  position 
of  the  towline  from  above.  Lastly,  the  heat  transfer  to  the  towline  was  modeled  under  a 
steady  state  position  over  time.  Due  to  the  dynamic  behavior  of  the  code,  and  since  no 
specific  situations  were  under  investigation,  the  general  behavior  of  the  towline  is 
represented. 

Although  not  shown  in  this  paper  (due  to  the  difficulty  of  showing  dynamic 
behavior  in  printed  material),  the  towline  can  also  be  modeled  as  a  pendulum  at  near  zero 
velocity  (zero  velocity  causes  a  slack  condition).  This  representation  reveals  a  damped 
oscillation  about  a  hanging  position,  with  curvature  in  the  line  as  the  body  gets  farther 
from  its  final  steady  state.  This  test,  which  was  used  to  test  the  MATLAB®  code  for 
reasonable  response,  is  an  easy  way  to  show  that  the  methods  used  work. 

4.1.  Initial  Conditions 

The  initial  conditions  for  the  towline  parameters  and  step  sizes  for  a  towline 
length  of  30  meters  were  set  as  follows: 
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ds  =5  m 
dt  =  .01  s 
dTL  =  .00127  in 

/u  =  .0096  k^~ 
m 

C  =1  1 
C  =04 
C  =1 

^D,B  1 

dB  =  .496  m 
Mb  =  1  kg 

where  most  of  these  values  were  taken  from  previous  work  by  Richardson.  The  value  for 
towline  mass  per  unit  length  (ju)  is  based  on  an  approximated  average  density  for  most 
steels  of -7600  kg/m  ,  and  is  a  function  of  both  density  and  diameter  (dn).  Schram  and 
Reyle  indicated  (Schram  and  Reyle,  1968:218)  that  the  value  of  ds/dt  should  be  around 
1000  to  get  good  results  and  ensure  stability  for  small  perturbations.  The  values  used 
here  seemed  to  be  stable  for  most  cases,  however.  Greater  changes  in  velocity  may 
require  different  values,  and  the  code  prints  an  error  statement  if  either  stability  cannot  be 
maintained  (Equations  3.20  and  3.21),  or  tension  becomes  negative. 

The  code  is  unable  to  accommodate  negative  tension  (i.e.,  compression)  in  the 
towline,  due  to  the  creation  of  imaginary  numbers  in  the  characteristic  directions 
(Equations  3.20  and  3.21).  This  is  consistent,  however,  with  the  assumptions  noted  by 
Crist  (see  Section  2.3)  in  the  fact  that  the  minimum  tension  a  line  hanging  freely  can  have 
is  zero.  Thus,  slack  conditions  take  place  at  zero  tension  (tension  cannot  equal  zero  in 
code  due  to  division  by  zero  errors,  so  slack  is  assumed  to  happen  at  near  zero  values). 

Unfortunately,  the  code  is  also  currently  unable  to  account  for  slack  conditions. 
This  is  because  the  characteristic  lines  become  zero  with  zero  tension.  A  line  of  slope 


64 


zero  will  never  be  able  to  transfer  data  between  timesteps  (see  Figure  3.3-1).  Thus, 
consideration  must  be  made  in  the  towed  body’s  drag  and  mass,  which  can  be  altered  to 
prevent  slack.  It  was  found  that  usually  either  reducing  the  mass,  or  increasing  the  towed 
body’s  drag  coefficient  or  frontal  area  (i.e.,  increase  the  drag)  prevented  this  condition. 
Although  Crist  notes  that  “the  effects  of  slack  cable  must  be  included  in  most  analyses” 
(Crist,  1970:73),  the  design  of  many  towed  bodies  are  made  to  prevent  slack  in  the  line, 
and  it  is  assumed  that  an  aircraft  in  combat  will  work  at  high  enough  speeds  such  that 
slack  becomes  nonexistent.  Line  slack  will  create  adverse  effects  such  as  line  jerk  that 
are  undesirable  anyway.  The  values  for  mass  and  drag  are  changed  in  some  of  the 
following  scenarios  to  prevent  this  effect,  and  any  future  work  should  note  the  restriction 
on  this  method. 

4.2.  Steady  State  Analysis 

Using  the  steady  state  values  from  Richardson’s  code,  the  towline  behavior 
should  remain  relatively  stable.  If  the  code  was  started  from  a  towline  hanging  position 
(either  back  or  down,  noting  an  initial  nonzero  angle  must  be  given  for  the  backward 
position),  the  behavior  should  eventually  reach,  or  maintain  (if  given  initial  steady  state 
values),  a  steady  state  value  close  to  that  from  Richardson’s  code. 

Indeed,  a  steady  state  value  is  found  that  is  consistent  with  that  from  Richardson. 
Using  the  initial  values  from  Richardson,  any  towline  movement  is  negligible,  and  a  final 
shape  can  be  seen  to  directly  overlap  the  one  from  his  code.  This  is  shown  in  a  50  m/s 
horizontal  velocity  simulation  plotted  in  Figure  4.2-1,  where  the  towline  was  started  at 
Richardson’s  steady  state  position. 
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Steady  State:  u  =  50m/s,  Cd  =  1 ,  m  =  1  kg 


Figure  4.2-1  Steady  State  Comparison 


It  should  be  noted  that  the  mesh  spacing  was  reduced  to  a  ds  value  of  1  meter  and 
a  dt  value  of  .001  seconds  to  increase  accuracy.  Note  the  very  slight  difference  between 
the  two  methods  near  the  towed  body.  This  difference  is  small  enough  to  be  negligible. 
Any  difference  is  most  likely  due  to  the  accuracy  of  the  method  of  characteristics,  which 
is  not  as  accurate  as  the  procedure  Richardson  used.  His  work  utilized  MATLAB®’s 
ode45  solver,  which  is  a  4th  order  Runge-Kutta  integration  scheme.  The  method  of 
characteristics  for  this  work  used  a  preset  mesh  with  straight  characteristic  lines,  which  is 
only  1st  order  accurate.  Thus,  some  small  discrepancies  may  be  found  between  the  two 
procedures.  As  mentioned,  however,  these  are  mostly  negligible. 

A  final  shape  can  be  seen  by  also  modeling  the  hanging  position.  As  mentioned 
before,  the  towline  was  started  at  a  constant  angle  of  (j)  =  0.01  radians  and  released.  The 
aircraft  was  kept  at  a  constant  speed  of  50  m/s,  and  the  simulation  took  place  over  50 
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seconds.  The  results  can  be  shown  in  Figure  4.2-2  through  Figure  4.2-6.  The  steady  state 
values  from  Richardson  are  included  on  these  plots  in  order  to  see  the  oscillations  about 
steady  state  and  model  the  damped  response.  The  towline  eventually  reaches  the  same 
steady  state  position  as  shown  previously  in  Figure  4.2-1.  Tension  is  included  in  these 
results  to  show  how  it  behaves  over  time.  Note  that  the  tension  is  always  less  at  the  body 
than  at  the  aircraft. 
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Steady  State:  u  =  50m/s,  Cd  =  1 ,  m  =  1kg,  time  =  Is 


Figure  4.2-3  Steady  State  Drop  t=ls 


Steady  State:  u  =  50m/s,  Cd  =  1 ,  m  =  1kg,  time  =  1  5s 


Figure  4.2-4  Steady  State  Drop  t=1.5s 
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Steady  State:  u  =  50m/s,  Cd  =  1 ,  m  =  1kg,  time  =  3s 


Figure  4.2-5  Steady  State  Drop  t=3s 


Steady  State:  u  =  50m/s,  Cd  =  1 ,  m  =  1kg,  time  =  50s 


Figure  4.2-6  Steady  State  Drop  t=50s 
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The  towline  initially  goes  down,  crossing  the  steady  state  value.  The  tension 
increases  until  this  value  starts  to  be  crossed.  At  this  point,  the  tension  decreases,  and  the 
towline  starts  to  slow  down.  What  happens  physically  during  this  time  is  that  the  drag 
forces  on  the  line  and  body  start  to  catch  up  to  the  gravitational  forces,  thus  reducing  the 
tension.  The  towline  then  goes  back  up  past  the  steady  state  value,  but  not  as  far.  This 
oscillation  happens  until  a  steady  state  value  is  reached.  This  value  is  very  close  to  that 
of  Richardson,  such  that  the  differences  are  negligible,  and  is  the  same  shape  as  that 
shown  in  Figure  4.2-1. 

4.3.  Single  Perturbation 

The  aircraft  was  set  at  an  initial  horizontal  velocity  of  100  m/s  in  a  similar  manner 
as  in  the  steady  state  analysis.  The  reference  frame  for  all  of  these  plots  is  in  the 
aircraft’s  initial  velocity  and  direction  (100  m/s  horizontal  velocity,  which  is  to  the  right 
in  the  x-y  plane).  A  way  to  visualize  this  is  to  think  of  being  in  a  chase  aircraft,  traveling 
at  a  constant  velocity  with  the  initial  velocity  of  the  aircraft  carrying  the  towline. 

The  code  is  set  such  that  the  acceleration  is  constant  over  a  predetermined  period 
of  time.  This  creates  a  jerk  effect  in  the  tension  (not  shown  here),  which  could  be 
reduced  by  ramping  the  acceleration.  Since  the  jerk  is  in  acceleration,  which  is  the 
second  derivative  of  position,  this  has  no  noticeable  effect  on  how  the  system  behaves  in 
space,  and  thus  is  mentioned  only  as  an  explanation  of  how  the  code  works. 

Horizontal  Perturbation. 

The  decoy  mass  was  set  to  4  kg  for  this  scenario.  The  first  perturbation  had  the 
aircraft  slow  down  to  70  m/s  over  a  period  of  1  second,  then  immediately  sped  back  up  to 
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y  (m)  y  (m) 


100  m/s  over  a  period  of  1  second.  The  horizontal  acceleration  is  shown  in  Figure  4.3-1 
through  Figure  4.3-7. 


-30  Horizontal  (1  s  accel,  1  s  start,  1  s  accel):  u  =  lOOm/s,  Cd  =  1 ,  m  =  4kg,  time  =  Os 
1  — i - 1 - 1 - 1 - 1 - 1 - r 
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Figure  4.3-1  Horizontal  Perturbation  t=0s 


-30  Horizontal  (Is  accel,  1  s  start,  Is  accel):  u  =  lOOm/s,  Cd  =  1 ,  m  =  4kg,  time  =  ,5s 
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-30  Horizontal  (Is  accel,  1  s  start,  Is  accel):  u  =  lOOm/s,  Cd  =  1 ,  m  =  4kg,  time  =  Is 


x  (m) 

Figure  4.3-3  Horizontal  Perturbation  t=ls 


-30  Horizontal  (Is  accel,  1  s  start,  Is  accel):  u  =  lOOm/s,  Cd  =  1 ,  m  =  4kg,  time  =  1.5s 
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-30  Horizontal  (Is  accel,  1  s  start,  Is  accel):  u  =  lOOm/s,  Cd  =  1 ,  m  =  4kg,  time  =  2s 


-30  Horizontal  (Is  accel,  1  s  start,  Is  accel):  u  =  lOOm/s,  Cd  =  1 ,  m  =  4kg,  time  =  3s 
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-30  Horizontal  (Is  accel,  1  s  start,  Is  accel):  u  =  lOOm/s,  Cd  =  1 ,  m  =  4kg,  time  =  20s 


As  can  be  seen,  the  entire  system  slows  down  and  the  body  drops  below  its  initial 
and  final  steady  state  value.  Upon  acceleration  back  to  100  m/s,  the  body  then  comes  up 
past  its  steady  state  value,  and  after  oscillating  about  this  value  settles  down  in  the 
original  steady  state  position. 

Vertical  Perturbation. 

This  scenario  also  used  a  4  kg  mass.  The  aircraft  was  accelerated  to  a  30  m/s 
vertical  velocity  over  1  second  of  acceleration.  This  velocity  was  kept  constant  until  2 
seconds,  at  which  point  it  was  accelerated  back  to  its  initial  vertical  velocity  of  0  m/s  over 
a  period  of  1  second.  The  results  can  be  seen  in  Figure  4.3-8  through  Figure  4.3-16. 
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+30  Vert  Movement  (1  s  accel,  2s  start,  Is  accel):  u  =  lOOm/s,  Cd  =  1 ,  m  =  4kg,  time  =  Os 
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+30  Vert  Movement  (Is  accel,  2s  start,  1  s  accel):  u  =  lOOm/s,  Cd  =  1 ,  m  =  4kg,  time  =  ,5s 


Figure  4.3-8  Vertical  Perturbation  t=Os 
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+30  Vert  Movement  (Is  accel,  2s  start,  Is  accel):  u  =  lOOm/s,  Cd  =  1 ,  m  =  4kg,  time  =  1.5s 


Figure  4.3-10  Vertical  Perturbation  t=1.5s 

+30  Vert  Movement  (Is  accel,  2s  start,  Is  accel):  u  =  lOOm/s,  Cd  =  1 ,  m  =  4kg,  time  =  2.5s 
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+30  Vert  Movement  (1  s  accel,  2s  start,  Is  accel):  u  =  lOOm/s,  Cd  =  1 ,  m  =  4kg,  time  =  3s 


Figure  4.3-12  Vertical  Perturbation  t=3s 

+30  Vert  Movement  (1  s  accel,  2s  start,  Is  accel):  u  =  lOOm/s,  Cd  =  1 ,  m  =  4kg,  time  =  3.5s 
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+30  Vert  Movement  (Is  accel,  2s  start,  Is  accel):  u  =  lOOm/s,  Cd  =  1 ,  m  =  4kg,  time  =  4s 
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+30  Vert  Movement  (Is  accel,  2s  start,  Is  accel):  u  =  lOOm/s,  Cd  =  1 ,  m  =  4kg,  time  =  4.5s 


Figure  4.3-14  Vertical  Perturbation  t=4s 
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+30  Vert  Movement  (1  s  accel,  2s  start,  Is  accel):  u  =  lOOm/s,  Cd  =  1 ,  m  =  4kg,  time  =  50s 


Figure  4.3-16  Vertical  Perturbation  t=50s 

Note  the  overshoot  as  the  system  returns  to  its  initial  vertical  velocity.  When 
applying  this  to  the  concern  of  the  towline  crossing  behind  the  engine,  assuming  that  the 
engine  plume  is  always  in  the  horizontal  direction,  the  towline  can  be  seen  to  cross  right 
down  the  middle  of  the  plume.  This  effect  is  not  desirable  when  concerned  about  heat 
transfer  to  the  line.  Since  the  towline  is  generally  not  deployed  exactly  behind  the 
engine,  and  behavior  represented  here  is  only  in  the  vertical  direction,  this  is  not 
necessarily  the  worst  case.  However,  it  is  one  that  should  be  avoided. 

Although  it  is  not  very  noticeable  in  these  plots,  the  towline  exhibits  a  curvature 
in  the  overshoot,  which  is  expected.  The  reason  for  the  small  amount  of  curvature  is  that 
the  mass  is  contributing  a  large  force  in  the  y-axes  due  to  weight.  The  inertia  from  the 
overshoot  is  what  causes  the  curvature,  but  since  gravity  is  always  acting  on  the  body,  the 
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force  in  the  v-axis  does  not  change  enough  to  create  a  significant  curving  effect. 
Essentially,  the  weight  serves  to  damp  the  curvature  from  inertia. 


Transverse  Perturbation. 


This  scenario  also  used  a  4  kg  mass.  The  aircraft  was  accelerated  to  a  30  m/s 
velocity  in  the  transverse  direction  (z-axis)  over  1  second  of  acceleration.  At  1  second, 
this  velocity  was  accelerated  back  to  its  initial  vertical  velocity  of  0  m/s  over  a  period  of 
1  second.  The  results  can  be  seen  in  Figure  4.3-17  through  Figure  4.3-26. 


+30  Transverse  (Is  accel,  1  s  start,  Is  accel):  u  =  50m/s,  Cd  =  1 ,  m  =  4kg,  t  =  Os 
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Figure  4.3-17  Transverse  Perturbation  t=0s 
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+30  Transverse  (Is  accel,  1  s  start,  Is  accel):  u  =  50m/s,  Cd  =  1 ,  m  =  4kg,  t  =  ,5s 


Figure  4.3-18  Transverse  Perturbation  t=0.5s 


+30  Transverse  (Is  accel,  1  s  start,  Is  accel):  u  =  50m/s,  Cd  =  1 ,  m  =  4kg,  t  =  Is 


Figure  4.3-19  Transverse  Perturbation  t=ls 
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+30  Transverse  (Is  accel,  1  s  start,  Is  accel):  u  =  50m/s,  Cd  =  1 ,  m  =  4kg,  t  =  1.5s 


Figure  4.3-20  Transverse  Perturbation  t=1.5s 


+30  Transverse  (Is  accel,  1  s  start,  Is  accel):  u  =  50m/s,  Cd  =  1 ,  m  =  4kg,  t  =  2s 


Figure  4.3-21  Transverse  Perturbation  t=2s 
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+30  Transverse  (Is  accel,  1  s  start,  Is  accel):  u  =  50m/s,  Cd  =  1 ,  m  =  4kg,  t  =  2.5s 


Figure  4.3-22  Transverse  Perturbation  t=2.5s 


+30  Transverse  (Is  accel,  1  s  start,  Is  accel):  u  =  50m/s,  Cd  =  1 ,  m  =  4kg,  t  =  3s 


Figure  4.3-23  Transverse  Perturbation  t=3s 
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+30  Transverse  (Is  accel,  1  s  start,  Is  accel):  u  =  50m/s,  Cd  =  1 ,  m  =  4kg,  t  =  3.5s 


Figure  4.3-24  Transverse  Perturbation  t=3.5s 


+30  Transverse  (Is  accel,  1  s  start,  Is  accel):  u  =  50m/s,  Cd  =  1 ,  m  =  4kg,  t  =  5s 


Figure  4.3-25  Transverse  Perturbation  t=5s 
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Figure  4.3-26  Transverse  Perturbation  t=20s 


Since  the  towline  was  not  initially  in  a  steady  state  position,  the  change  in  both 
the  x-y  and  the  z-y  planes  can  be  seen.  Although  not  explicitly  clear  here,  the  changes  in 
the  transverse  direction  affect  the  changes  in  the  vertical  direction  and  vice  versa.  The 
curvature  of  the  line  is  more  distinct  in  these  perturbations,  since  the  weight  only  works 
in  thc  v-axis.  Thus,  the  inertia  in  the  z-axis  becomes  significant,  causing  curvature.  Of 
interesting  note  is  that  although  the  line  goes  past  its  steady  state  value,  it  does  not  come 
back  and  cross  this  position,  which  is  most  likely  due  to  the  small  amount  of  overshoot. 
Since  the  force  from  weight  is  what  brings  the  line  back  to  steady  state,  and  this  force  is 
relatively  small  due  to  the  small  angle  the  towline  goes  past  its  steady  state  value,  the 
towline  approaches  its  final  position  at  a  slow  rate,  thus  producing  little  to  no  overshoot 
when  it  comes  back.  A  greater  acceleration  in  the  transverse  direction  would  produce  a 
noticeable  overshoot  on  the  return  here. 
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4.4.  Multiaxial  Perturbation 


Using  the  same  values  for  the  previous  simulations  with  a  4  kg  mass,  the  last  two 
maneuvers  from  Section  4.3  were  combined  to  produce  an  upward  and  rightward 
movement.  Here,  however,  the  transverse  velocity  did  not  begin  to  accelerate  back  to 
zero  until  the  2  second  mark.  The  results  can  be  seen  in  Figure  4.4-1  through  Figure 
4.4-10,  noting  the  scaling  in  the  z-y  plane,  which  changes  so  that  the  results  can  be  more 
easily  seen. 


+30  Vert  &  +30  Transverse  Movement  (Is  accel,  2s  start,  Is  accel):  u  =  lOOm/s,  Cd  =  1 ,  m  =  4kg,  time  =  Os 
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Figure  4.4-1  Multiaxial  Perturbation  t=0s 
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+30  Vert  &  +30  Transverse  Movement  (Is  accel,  2s  start,  Is  accel):  u  =  100m/st  Cd  =  1 ,  m  =  4kg,  time  =  Is 


Figure  4.4-2  Multiaxial  Perturbation  t=ls 
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+30  Vert  &  +30  Transverse  Movement  (Is  accel,  2s  start,  Is  accel):  u  =  100m/st  Cd  =  1 ,  m  =  4kg,  time  =  3s 
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Figure  4.4-4  Multiaxial  Perturbation  t=3s 
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+30  Vert  &  +30  Transverse  Movement  (Is  accel,  2s  start,  Is  accel):  u  =  lOOm/s,  Cd  =  1 ,  m  =  4kg,  time  =  3.5s 


Figure  4.4-6  Multiaxial  Perturbation  t=3.5s 
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+30  Vert  &  +30  Transverse  Movement  (Is  accel,  2s  start,  Is  accel):  u  =  100m/st  Cd  =  1 ,  m  =  4kg,  time  =  4s 
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Figure  4.4-8  Multiaxial  Perturbation  t=4s 


+30  Vert  &  +30  Transverse  Movement  (Is  accel,  2s  start,  Is  accel):  u  =  lOOm/s,  Cd  =  1 ,  m  =  4kg,  time  =  5s 
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Figure  4.4-9  Multiaxial  Perturbation  t=5s 
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The  towline  overshoots  when  the  aircraft  returns  to  its  steady  state  value  (note  the 
aircraft’s  final  position  at  z  =  60m  and  y  =  60m  with  the  body  overshooting  these  values). 
It  appears  to  cross  right  behind  the  aircraft.  Applying  this  to  the  concern  of  heat  transfer 
to  the  towline,  this  case  is  potentially  worse  than  the  purely  vertical  perturbation.  The 
towline  crosses  behind  the  aircraft  in  both  the  y  and  z  axes,  which  could  put  it  right  in  the 
engine  plume  if  the  engine  is  off  center  from  the  towline  attachment. 

4.5.  Heat  Transfer 

All  the  previous  analyses  were  assumed  to  be  in  standard  atmosphere  and 
pressure,  thus  keeping  a  constant  air  density  with  constant  line  temperature.  To  analyze 
the  change  in  towline  temperature  and  air  density  about  the  towline,  a  routine  was  written 
( Tempset.m )  to  calculate  these  values  at  every  point  in  space  and  time  and  is  included  in 
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Appendix  C.  The  code  was  written  with  the  intention  of  given  data  being  placed  into  an 
ambient  temperature  array  as  global  or  otherwise  referenced  values. 

In  the  analysis  here,  however,  it  was  assumed  that  the  ambient  temperature  was 
held  constant  at  1000  K,  and  although  the  array  is  a  function  of  time  and  three 
dimensional  space  (relative  to  the  initial  positions  of  the  system),  all  of  its  values  are  the 
same.  This  can  be  easily  changed  by  adding  known  data.  Both  the  heat  transfer  to  the 
towline,  and  the  air  density  about  the  towline,  were  analyzed  under  this  constant  ambient 
temperature. 

The  vertical  velocity  perturbation  in  Section  4.3  was  run  in  the  same  manner  as 
before  over  a  period  of  150  seconds.  A  1  kg  mass  was  used  for  this  analysis,  however, 
since  a  lower  mass  will  produce  more  curvature  in  the  towline  (body  weight  works  to 
damp  any  curvature).  The  vertical  maneuver  was  used  in  order  to  model  how  heat 
transfer  is  affected  by  a  change  in  perpendicular  velocity  over  the  towline.  The  aircraft 
was  started  at  a  100  m/s  horizontal  velocity  and  was  given  an  initial  vertical  acceleration 
over  a  period  of  1  second  to  a  final  vertical  velocity  of  30  m/s.  This  was  held  for  1 
second.  At  the  2  second  mark,  the  same  procedure  is  done  in  reverse  to  slow  the  aircraft 
down  to  0  m/s  vertical  velocity.  As  mentioned,  however,  this  time  it  was  placed  in  an 
ambient  temperature  field  of  1000  K.  The  towline  temperature  was  given  an  initial 
temperature  of  288.2  K  (value  for  standard  atmosphere  and  pressure). 

The  first  plot  (Figure  4.5-1)  shows  the  temperature  in  the  towline  as  a  function  of 

time. 
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Temp  Change  Over  Time  (Tamb  =  1000  K) 


Figure  4.5-1  Temp  Change  Over  Time 


Note  that  the  line  temperature  asymptotically  approaches  the  ambient  temperature  value 
of  1000  K.  This  creates  a  near  linear  change  initially  and  a  much  more  gradual  change  as 
time  goes  on.  The  reason  that  the  line  temperature  at  the  aircraft  increases  faster  is  due  to 
the  towline  shape  in  its  steady  state  value.  This  is  shown  in  Figure  4.5-2. 
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Towline  Shape  at  Steady  State  (t  =  150s) 


Figure  4.5-2  Towline  Shape  at  Final  Steady  State 

Note  the  shape  of  the  towline,  which  allows  for  a  greater  perpendicular  velocity  over  the 
line  near  the  aircraft  (right  side  of  plot)  versus  near  the  body  (left  side  of  plot)  since  the 
system  is  traveling  at  100  m/s  to  the  right. 

Also  of  note  is  the  increased  droop  in  the  towline.  The  final  droop  is  about  0.93 
meters,  whereas  the  initial  droop  is  about  0.64  meters.  The  initial  position  of  the  towline 
before  the  maneuver  is  shown  in  Figure  4.5-3. 
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Towline  Shape  at  Tamb  =  288. 2K 


Figure  4.5-3  Towline  Shape  at  Initial  Steady  State 
This  increased  droop  is  due  to  the  reduction  in  drag  forces  as  a  result  of  reduced  air 
density.  The  air  density  is  shown  in  Figure  4.5-4. 
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Air  Density  Change  Over  Time  (TQmb  =  1000K) 


Note  how  the  density  changes  as  a  function  of  temperature.  A  greater  temperature  causes 
a  lower  air  density  under  constant  pressure. 

To  illustrate  the  initial  rate  of  heat  transfer,  the  towline  temperature  is  plotted  for 
the  first  5  seconds  in  Figure  4.5-5.  The  towline  shape  can  be  seen  over  this  time  period  in 
Figure  4.3-8  through  Figure  4.3-15.  Note  that,  over  this  time  period,  the  towline  is 
traveling  faster  near  the  body  than  near  the  aircraft.  Thus,  there  is  greater  heat  transfer 
near  the  body  than  near  the  aircraft.  The  two  values  cross  around  10  seconds,  due  to  a 
slightly  different  slope  in  the  steady  state  position.  It  should  be  noted  that  the  slope  is  the 
heat  transfer  rate. 
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lamp  Change  On  Tne  "  1000  K) 


The  change  in  the  heat  transfer  rate  at  the  aircraft  is  noticeable  at  the  three 
inflection  points  for  acceleration  ( 1  second,  2  seconds,  and  3  seconds).  The  first  second 
has  the  aircraft  increasing  in  velocity  (curved  heat  transfer  slope).  The  second  second  has 
the  aircraft  at  a  constant  velocity,  but  higher  than  its  initial  and  final  velocities  (linear 
heat  transfer).  The  third  second  has  the  aircraft  slowing  down  to  its  final  velocity  (curved 
heat  transfer).  After  this  final  acceleration  is  complete  and  the  system  approaches  steady 
state  again,  the  heat  transfer  near  the  aircraft  settles  to  an  essentially  linear  behavior.  The 
curvature  between  3  and  4  seconds  is  due  to  the  system  approaching  its  steady  state.  The 
heat  transfer  at  the  body  becomes  a  little  more  complex  due  to  its  oscillations  about  the 
steady  state  value,  but  follows  expected  behavior.  Note  the  slightly  higher  heat  transfer 
rate  at  the  body  on  the  right  side  of  the  plot. 
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V:  Conclusions  and  Recommendations 


This  paper  utilized  a  variety  of  research  to  develop  a  reliable  code  to  model 
towline  behavior  under  transient  conditions.  When  compared  to  past  work,  it  has  been 
shown  to  accurately  predict  the  position  of  the  towline  under  any  transient  maneuvers. 
This  code  can  also  accurately  model  the  heat  transfer  to  the  line  based  on  given  initial 
data  within  the  assumptions  given  (i.e.,  not  heat  transfer  due  to  parallel  airflow).  The 
theory  is  presented  clearly  and  the  code  works  effectively. 

5.1.  Conclusions 

A  few  conclusions  can  be  drawn  from  the  results  shown  here.  The  first 
conclusion  is  that  the  method  of  characteristics  gives  a  good  approximation  for  the 
towline  behavior  in  both  transient  and  steady  state  analysis.  This  was  shown  in 
comparison  with  past  work,  as  well  as  the  modeling  of  expected  behaviors  such  as  the 
pendulum  test. 

The  second  conclusion  is  that  this  method  cannot  analyze  slack  conditions.  Since 
most  towed  decoys,  which  are  the  focus  of  this  research,  are  designed  to  prevent  slack, 
this  method  should  be  able  to  analyze  most  cases  of  concern.  Although  the  slack 
conditions  cannot  be  analyzed,  the  towline  acts  as  it  should  when  compared  to  previous 
work.  As  a  result,  special  consideration  needs  to  be  made  regarding  the  towed  body’s 
drag  and  mass,  which  must  be  adjusted  in  order  to  prevent  slack. 

The  third  conclusion  regards  the  coupling  of  motion.  Previous  research  has 
indicated  that  disturbances  in  one  plane  generate  disturbances  in  another  (Schram  and 
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Reyle,  1968:219).  This  is  seen  in  the  results  here  as  well  except  for  the  case  of  a  straight 
towline  undergoing  maneuvers  purely  in  the  plane  in  which  it  starts  (i.e.,  a  purely  vertical 
maneuver  from  straight,  steady  state  flight).  Previous  work  indicated  this  exception  as 
well  (Schram  and  Reyle,  1968:219),  so  the  effect  is  not  unexpected.  This  coupling  of 
motion  is  due  to  the  change  in  drag  being  put  on  the  line  and  body  due  to  the  change  in 
total  velocity,  which  affects  the  towline’s  behavior  in  multiple  axes. 

The  fourth  conclusion  is  that  the  heat  transfer  to  the  towline  appears  linear  over 
short  periods  of  time.  A  gradual  leveling  off  over  longer  periods  of  time  happens,  but  for 
shorter  periods,  the  change  appears  linear.  This  is  expected  for  large  differences  in 
temperature,  since  Newton’s  law  of  cooling  says  that  the  rate  of  heat  transfer  is 
dependent  on  the  difference  in  temperatures.  Thus,  a  smaller  difference  produces  a 
smaller  change. 

5.2.  Future  Work 

Future  work  entails  the  development  of  a  better  way  to  approximate  the  Nusselt 
number,  analyzing  the  boundary  layer  on  the  towline,  modeling  the  effects  of  wake  from 
the  aircraft,  modeling  slack  conditions  in  the  line,  the  use  of  different  mesh  schemes, 
analyzing  the  towed  body’s  forces,  and  comparing  to  real  results  for  specific  applications 
(i.e.,  specific  maneuvers).  In  addition,  in  order  to  analyze  heat  transfer,  a  temperature 
field  as  a  function  of  time  and  space  must  be  referenced  from  known  data. 

As  alluded  to  in  previous  sections,  the  method  derived  is  limited  in  its  ability  to 
model  heat  transfer  to  the  towline  due  to  the  Nusselt  number.  Future  research  should 
entail  a  manner  in  which  to  find  the  Nusselt  number  at  differing  geometries,  such  that  the 
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heat  transfer  to  the  line  can  be  approximated  for  both  parallel  and  perpendicular  flow. 
This  may  require  boundary  layer  analysis  on  the  towline. 

The  drag  and  heat  transfer  from  airflow  parallel  to  the  towline  were  both  assumed 
to  be  negligible  for  this  work.  One  of  the  primary  causes  of  this  assumption  is  the 
formation  of  a  boundary  layer,  which  increases  in  thickness  along  the  towline.  There  is 
no  sure  way  to  analyze  the  boundary  layer,  however,  because  flow  separation  will  start  to 
occur  along  different  points  of  the  line.  Even  though  flow  separation  occurs,  it  was 
assumed  still  that  the  length  size  of  the  turbulence  was  much  greater  than  that  of  each 
incremental  length  of  the  towline.  Due  to  this,  all  the  airflow  becomes,  essentially, 
perpendicular,  and  thus  a  perpendicular  airflow  analysis  may  be  sufficient.  More  work 
should  be  done,  however,  to  model  the  formation  of  the  boundary  layer,  its  effects  on 
drag  and  heat  transfer,  and  the  result  of  turbulent  effects  due  to  flow  separation. 

All  flow  was  considered  to  be  laminar  at  every  local  point  along  the  towline.  This 
is  not  necessarily  true,  especially  since  this  system  is  being  towed  behind  an  aircraft.  An 
analysis  of  the  turbulent  wake  from  the  aircraft  would  determine  if  this  would  affect  any 
of  the  data. 

Slack  conditions  cannot  be  analyzed  with  this  model.  That  is  because  the 
characteristics  used  to  step  data  from  one  timestep  to  the  next  are  functions  of  tension.  A 
zero  tension  value  prevents  the  transfer  of  data  by  this  method.  Thus,  another  method,  or 
some  alteration  of  the  one  used  here,  is  necessary  to  allow  for  slack  in  the  line.  This 
method  could  be  added  to  the  current  code  as  an  if  statement,  called  up  only  during  slack 
conditions. 
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Although  the  first  order  scheme  of  Courant,  Isaacson,  and  Rees  (Ames,  1965:447) 
used  here  for  the  method  of  characteristics  provides  a  very  simple  mesh  setup,  their’ s  is 
only  first  order  accurate.  A  better  setup  for  the  method  of  characteristics  involves  a 
second  order  solution  based  on  Hartree’s  hybrid  method,  as  outlined  by  Ames  (Ames, 
1965:445-448).  The  technique  used  in  this  paper  is  essentially  an  approximation  of  this 
hybrid  method.  Although  probably  not  necessary,  the  alternate  method  Ames  lays  out 
would  provide  a  more  accurate  approximation,  especially  at  larger  step  sizes.  A  larger 
step  size  would  reduce  computing  time.  Unfortunately,  it  also  requires  iterative 
procedures  for  solutions,  which  increases  computing  time.  The  payoffs  or  advantages  of 
using  this  scheme  are  unknown.  However,  if  computing  time  becomes  an  issue,  more 
should  be  looked  into  using  the  second  order  scheme  instead.  Also,  as  was  shown  in  the 
steady  state  analysis  of  Section  4.2,  there  is  a  slight  discrepancy  between  the  1st  order 
method  used  here,  and  the  4th  order  method  used  by  Richardson.  Using  this  second  order 
method  from  Ames  should  reduce  the  discrepancy.  These  errors  are  small,  however,  and 
the  work  necessary  to  rewrite  all  the  code  may  not  be  worth  the  slightly  increased 
accuracy. 

One  assumption  that  was  made  involved  the  orientation  of  the  body  being  always 
facing  directly  into  the  freestream.  This  is  not  a  bad  assumption,  but  future  work  should 
involve  a  better  modeling  of  the  forces  on  the  body.  Depending  on  the  body’s  weight  and 
drag,  it  can  have  a  very  significant  effect  on  the  towline  behavior.  Since  the  body  shape 
was  unknown  and  only  assumptions  were  made  for  this  work,  future  analyses  using 
actual  body  properties  will  help  produce  very  useful  data. 


101 


As  an  all  theoretical  work,  data  from  this  code  should  be  compared  to  real  life 
data.  There  are  many  real  life  simulations  available.  A  lot  of  these  (including  other 
methods  used  to  analyze  this  type  of  problem)  can  be  found  referenced  in  Kang  and 
Latorre  (Kang  and  Latorre,  1991 :5-6).  Future  work  should  entail  comparing  the  results 
from  this  code  to  specific  scenarios  from  experiment  in  order  to  predict  behavior  under 
the  known  conditions. 

Other  possible  work  involves  some  small  additions  to  the  code,  both  of  which 
would  be  quite  simple.  These  include  utilizing  a  ramp  function  for  acceleration  to  reduce 
the  “jerk”  in  the  tension,  and  adding  more  possible  maneuvers  for  each  simulation 
(currently  only  2  are  allowed). 

These  conclusions  and  recommendations  can  be  further  substantiated  with  more 
research.  For  now,  however,  the  attached  code  is  sufficient  for  modeling  transient 
maneuvers  in  aircraft  behavior,  as  well  as  heat  transfer  to  the  towline  in  a  known 
temperature  field. 
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Appendix  A:  Development  of  Angular  Acceleration  Term 


Schram  (Schram,  1968:  Appendix  B)  developed  a  manner  in  which  to  calculate 
the  three-dimensional  rotational  equations  of  motion  of  the  body.  A  simplified  version  of 
a  manner  in  which  to  find  angular  acceleration  can  be  based  off  of  Equation  3.27,  where 
it  is  noted  that  the  angular  terms  become 

os^(A-A)] 

for  the  lower  boundary.  This  can  be  added  to  the  tension  at  the  body  to  account  for 
angular  acceleration  at  the  body. 
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Appendix  B:  MATLAB®  Code  -  Primary  Code 


%Theoretical  Modeling  of  the  Transient  Effects  of  a  Towline 
%Using  the  method  of  characteristics  (June  2006) 

%Created  by  Ensign  Christopher  A.  Hill,  USN 

%Under  the  advice  of  Dr.  Ralph  Anthenien 

%Master  Thesis  work.  Air  Force  Institute  of  Technology 

%For  more  information  contact  Christopher  Hill  at  612-532-6068 

%or  Ralph  Anthenien  at  937-255-3636  x4643 

%Primary  Code 

%all  arrays  of  form  f (position, time) 

clear  all 
clc 

%Amatrix  is  a  function  of  the  form  Amatrix (phi , theta) ,  which  forms  the 
%transf ormation  matrix  [A]  to  convert  from  space  to  towline  coordinates 
and 

%vice  versa. 

set  =  set  values;  %sets  variables  for  maneuvers 
params  =  line  params;  %brings  in  initial  conditions 

%Set  all  initial  values: 

t=0;  %sets  initial  time 

ds  =  set.ds;  %change  in  length  (meters) 
dt  =  set.dt;  %change  in  time  (sec) 
time  =  set. time;  %sec  (total  time) 

L  =  params. LL; 

mu  =  params. mu;  %mass  of  towline  per  meter 
g  =  params.g; 

Wt  =  params. Wt;  %Weight  of  towline  per  meter  [N] 

Mb  =  params.mD; 

WB  =  params.WB; 

Vol  =  set . ds*pi*params . dLA2 /4  ;  %total  volume  of  towline 
Area  =  set . ds*pi*params . dL;  %total  surface  area  of  towline 
C  =  params. Vx;  %m/s  -  constants  (velocity)  of  a/c 
D  =  params. Vy; 

E  =  params.Vz; 

upert  =  set.upert;  %perturbations  in  velocity  of  a/c 
vpert  =  set.vpert; 

wpert  =  set.wpert;  %limits  in  these  values 

tpert  =  set.tpert;  ^constant  acceleration  until  reaches  final 
perturbation  value  at  tpert 

tstart2  =  set.tstart2;  %time  to  start  second  perturbation 

upert2  =  set.upert2; 

vpert2  =  set.vpert2; 

wpert2  =  set.wpert2; 

tpert2  =  set.tpert2; 
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can  be 


%Set  values  for  finding  heat  transfer: 

Cp  =  450;  %approx  value  for  most  steels  at  300K  - 
changed  to  a  fnct  of  temp 

Pr  =  0.7;  %Prandtl  number 
Temp_set  =  set . CT*ones ( 15 , 15 , 2 , set . time/set . dt+1 ) ;  %Ambient 
temp  field:  Temp  set(x,y,z,t)  -  currently  uses  const,  temp  value 
'set.CT'  -  value  of  1  is  initial  position/time 

GS  =  .1;  %Spacial  grid  size  spacing  for  Temp  set  is  1/GS 
meters.  Thus  GS  =  10  means  grid  size  for  Temp  set  is  .1  meters 

%Set  the  offset  from  the  a/c  at  initial  conditions  for  where 

the 

%values  begin  for  Temp  set  matrix  (+  is  up  and  -  is  down) : 
xoffset  =  -30; 
yoffset  =  -30; 
zoffset  =  0; 

%  Setup  Initial  Conditions  (this  is  to  find  steady  state) 

%  +1  is  to  account  for  the  zero  point  on  the  line  (i.e.  a  '1' 
position/ time 

%  refers  to  a  position/time  of  0,  etc.) 

%  note:  we'll  say  that  t=l  is  initial  given  values,  and  t=2  is  upper 
limit 

%  conditions 

%  note:  everything  is  oriented  about  the  b  position  in  fig.  3.3-1  of 
paper 

%  ( ' i '  position  in  code  refers  to  b  position  in  paper) 

%Initialize  arrays  for  data: 

Phi  =  set . Phi*ones (L/ds+1 , time/dt+1 ) ;  %set  an  initial  value  for  Phi 
for  tO  -  affects  possible  negative  tension  at  initial  values 
Theta  =  set . Theta*ones (L/ds+1 , time/dt+1 ) ; 

U  =  zeros (L/ds+1 , time/dt+1 ) ; 

W  =  zeros (L/ds+1 , time/dt+1 ) ; 

V  =  zeros (L/ds+1, time/dt+1) ; 

FX  =  zeros (L/ds+1 , time/dt+1 ) ;  %note:  per  meter! 

FY  =  zeros (L/ds+1 , time/dt+1 ) ; 

FZ  =  zeros (L/ds  +  1 , time/dt+1 )  ; 

T  =  zeros (L/ds+1 , time/dt+1 ) ; 
u  =  zeros (L/ds+1 , time/dt+1 ) ; 
v  =  zeros (L/ds+1 , time/dt+1 ) ; 
w  =  zeros (L/ds+1, time/dt+1) ; 
x  =  zeros (L/ds+1, time/dt+1) ; 
y  =  zeros (L/ds+1 , time/dt+1 ) ; 
z  =  zeros (L/ds+1, time/dt+1) ; 

%These  two  arrays  are  all  set  to  the  same  value,  with  the 
understanding 

%that  the  non  t=l  array  positions  will  be  changed  later 
Temp  =  set . LTemp*ones (L/ds+1 , time/dt+1 ) ;  %sets  entire  line  temp  to 
be  initial  line  temp 

rho  =  357 . 88* ( ( set . LTemp+set . CT) /2 ) A ( - 
1 . 0041) *ones (L/ds+1, time/dt+1) ;  %sets  air  density  along  line  based  on 
initial  air  and  line  temps 

%Set  velocities  down  entire  line  =  constant  vel.  values  of  a/c 
(spatial) 
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u  ( : ,  1 ) 

II 

a 

v ( : , 1) 

ii 

u 

w  (  :  , 1) 

w 

ii 

%The  code  between  the  lines  takes  in  previous  steady  state  data. 

%One  can  comment  out  this  code  if  they  would  rather  start  from 
%non-steady  state  cases. 

O, _ _ _  _ _ _ _ _  _  _  _  _ _ _ _ _  _  o 

o  o 

%%The  following  sets  up  initial  conditions  by  interpolation  from 
%%previous  tension  data  (since  we  have  more  data  pts  here) .  It  should 
%%be  noted  that  the  previous  data  plots  the  towed  body  at  the  zero 
position 

%%and  follows  the  towline  shape  in  the  positive  x-y  direction  up  to  the 
%%aircraft.  We  need  to  shift  everything  down  and  to  the  left  so  that 
the 

%%aircraft  is  at  zero  and  the  towline  hangs  down  and  to  the  left. 

%  Create  matrices  to  store  values: 
dx  =  zeros (L/ds+1 , time/dt+1 ) ; 
dy  =  zeros (L/ds+1 , time/dt+1 ) ; 
dz  =  zeros (L/ds+1 , time/dt+1 ) ; 

%  Find  steady  state  values: 

[yl  y2 ]  =  towlinecomp; 

x  init  =  y2 ( : , 1 ) ;  %spatial  x  pos  of  steady  state  line  wrt  arc 

length 

y_init  =  -y2(:,2);  %spatial  y  pos  of  ss  line  wrt  arc  length... 
negative  due  to  sine  convention 

z  init  =  y2(:,3);  %spatial  z  pos  of  ss  line  wrt  arc  length 
dx  init  =  -y2(:,4);  %spatial  dx/ds  of  ss  line...  negative  due 
to  sine  convention 

dy_init  =  -y2(:,5);  %spatial  dy/ds  of  ss  line...  negative  due 
to  sine  convention 

dz_init  =  y2(:,6);  %spatial  dz/ds  of  ss  line 

T  init  =  y2 ( : , 7 ) *params . TO ;  %tension  in  ss  line...  originally 
non-dim,  thus  multiplied  by  TO 

Last  =  length (x  init);  %references  last  pt  in  initial  data  (same  for 
all  initial  matrices) 

e  =  (Last-1) / (length (x (:, 1) ) -1) ;  %ratio  for  interpolation  (same  for  all 
matrices)  ...  -1  used  b/c  1st  point  is  really  at  0  (no  0  index) 

%The  following  puts  the  initial  values  into  the  initial  matrices 
(swapped 

%for  the  convention  in  this  program) : 
for  i=l:l:L/ds+l 

%This  sets  up  the  indices  for  conversion  from  previous  coordinates 

for 

%this  program  using  interpolation  (A  is  previous  point,  B  is  next 
point 

%... these  values  are  subtracted  from  the  'Last'  position  in  order 
to 

%reverse  the  index) : 

A  =  round ( (i-1) *e) ; 
if  (i-1) *e  >  A 
B  =  A+l ; 
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elseif  A  ==  0 
B  =  A+l; 

else 

B  =  A; 

A  =  A-l; 

end 

%These  are  all  in  spatial  coordinates: 

x(i,l)  =  [ (x  init (Last-B) * ( ( i — 1 ) *e-A)  +  x  init (Last-A) * (B- (i- 
l)*e))*L  -  x  init ( length (x  init) ) *L] ; 

y(i,l)  =  -[ (y_init (Last-B) *( (i-1) *e-A)  +  y_init (Last-A) * (B- (i- 
1) *e) ) *L  -  y_init (length (y_init) ) *L]  ; 

z(i,l)  =  -[ (z_init (Last-B) *( (i-1) *e-A)  +  z_init (Last-A) * (B- ( i- 
1) *e) ) *L  -  z  init ( length ( z  init) ) *L] ; 

dx(i,l)  =  dx  init (Last-B) *( (i-1) *e-A)  +  dx  init (Last-A) * (B- (i- 

1) *e)  ; 

dy(i,l)  =  dy_init (Last-B) *( (i-1) *e-A)  +  dy_init (Last-A) * (B- (i- 

1) *e)  ; 

dz(i,l)  =  dz  init (Last-B) *( (i-1) *e-A)  +  dz  init (Last-A) * (B- (i- 

1) *e)  ; 

%This  finds  the  initial  angles: 

Phi(i,l)  =  atan (dy ( i , 1 ) /dx ( i , 1 ) ) ; 

Theta(i,l)  =  -asin (dz (i, 1 ) ) ; 


%  for  i=0:l:L/ds  %comment  out  this  line  if  using  above  for  loop, 
otherwise 
%  uncomment 

O, 

o 

%  %This  sets  up  the  values  for  a  non-steady  state  case,  setting  the 

%  %towline  positions  to  be  directly  backward  from  the  aircraft. 

These 


%  %lines  should  be  commented  out  if  using  the  steady  state  values 

above : 


%  if  i  ~=  0 

%  x (i, 1) 

%  y  (i,  1) 

%  z (i, 1) 

%  end 

o, _ 

o 


x(i,l)  -  ds*cos (Theta (i, 1) ) *cos (Phi (i, 1) ) ; 
y  ( i , 1 )  -  ds*cos (Theta (i, 1) ) *sin (Phi  (i, 1) ) ; 
z(i,l)  -  ds*sin (Theta (i, 1) ) ; 


O, 

o 


%This  converts  velocity  to  towline  coordinates: 

Spatial_vel  =  [u(i,l);  v(i,l);  w(i,l)]; 

Tow  vel  =  Amatrix (Phi (i, 1) , Theta (i, 1) ) *Spatial  vel; 

U(i,l)  =  Tow_vel(l); 

V(i,l)  =  Tow_vel(2); 

W(i,l)  =  Tow_vel(3); 

%The  following  is  used  to  find  heat  transfer  and  density  in  ambient 
temp  at  next  timestep. 

%Due  to  computing  power,  it  can  be  commented  out  when  not  analyzing 
heat  transfer. 

Tempset; 

if  i  ~=  L/ds 
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%This  sets  up  the  aerodynamic  forces  on  the  line 
F  =  Towf orce (U ( i , 1 ) , V ( i , 1 ) , W ( i , 1 ) , rho ( i , 1 ) ) ; 

FX(i,l)  =  -F(l);  %note:  per  meter 

FY(i,l)  =  -F(2);  %notation  is  in  positive  X,Y,Z  axes 
FZ  ( i , 1 )  =  -F (3) ; 

end 

end 


%Find  drag  forces  on  the  body  and  tension  at  lower  end: 

%No  acceleration  terms  are  included  since  it  is  assumed  this  is 
%starting  from  rest. 

F  = 

Bodyforce (u(L/ds+l,l) ,v(L/ds+l,l) ,w(L/ds+l,l) , rho (L/ds+1, 1) ) ;  %Set 
forces  on  body  in  spatial  axes 
Fx  =  - F ( 1 ) ; 

Fy  =  -F (2) -WB; 

Fz  =  -F  (3) ; 

FTL  = 

Amatrix ( Phi (L/ds+1 , 1 ) , Theta (L/ds+1 , 1 ) ) * [Fx; Fy; Fz ] ; %Convert  to  towline 
coordinates 

T(L/ds+l,l)  =  -FTL (2 ) ; %Tension  is  directed  only  along  the 
towline  Y-axis 

%Set  aerodynamic  forces  on  line  at  lower  end: 

F  =  Amatrix (Phi (L/ds+1, 1) , Theta (L/ds+1, 1) )* [F (1)  ;  F (2 );  F (3)  ]  ; 
%Change  to  towline  coordinates 

FX (L/ds+1, 1)  =  -F(l);  %negative  sign  due  to  drag  being  in  same 
direction  as  velocity 

FY (L/ds+1, 1)  =  -F (2) ; 

FZ (L/ds+1, 1)  =  -F (3) ; 

%Set  tension  one  step  up  line 

T (L/ds, 1)  =  T (L/ds+1, 1) +ds* (- 
FY (L/ds,  1) +Wt*sin (Phi (L/ds, 1) ) *cos (Theta  (L/ds, 1) ) )  ; 

%Set  tension  up  rest  of  line 
for  i=L/ds-l : -1 : 1 

T  ( i , 1 )  =  T (i+2, 1) +2*ds* (- 
FY (i+1, 1) +Wt*sin (Phi (i+1, 1) ) *cos (Theta (i+1, 1) ) ) ; 
end 

for  t=l : 1 : time/dt;  %t=l  represents  tO  (1st  timestep)  -  we're  finding 
values  for  next  timestep  based  on  this  timestep  value 

%This  sets  up  the  perturbation  velocities  over  time  based  on 
constant 

%acceleration . . .  t  is  used  because  we're  concerned  with  the 
behavior 

%at  time  t+1  (i.e.,  assume  all  perturbation  velocities  are  0  at 

tO)  , 

%thus  the  total  change  in  time  up  to  that  point  will  be  t*dt: 

%First  perturbation  (u/v/wstar  values  set  twice  since  u/v/wstar 
%values  are  changed  in  the  second  perturbation,  but  are  based 
on 
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%constant  values  of  u/v/wstar  -  the  ' 1 '  values  -  since  the 

second 

%perturbation  happens  after  the  first  one) : 
if  t*dt  <=  tpert 

ustarl  =  upert*t*dt/tpert; 
vstarl  =  vpert*t*dt/tpert; 
wstarl  =  wpert*t*dt/tpert; 
ustar  =  ustarl; 
vstar  =  vstarl; 
wstar  =  wstarl; 

end 

%Second  perturbation: 

if  t*dt  >=  tstart2  &  t*dt  <=  tstart2+tpert2 

tO  =  t-tstart2/dt;  %timestep  that  this  starts  at 
ustar  =  ustarl  +  t0*upert2*dt/tpert2 ; 
vstar  =  vstarl  +  t0*vpert2 *dt/tpert2 ; 
wstar  =  wstarl  +  t0*wpert2 *dt/tpert2 ; 

end 

for  i=l:l:L/ds+l 

%Characteristic  directions: 

Fa  =  sqrt (T (i, t) /mu)  ; 

Fb  =  -sqrt (T (i, t) /mu) ; 
if  i==L/ds+l 

Fb=-sqrt (T (i, t) /Mb) ;  %sets  different  characteristic 
value  for  lower  bdry 
end 

%This  checks  for  stability  of  mesh: 
if  (abs (Fa) *dt/ds)  >=  1 

fprintf('Out  of  boundary ! \nTry  changing  ds  &  dt  (set 
ds/dt  to  at  least  >=  1000) . \nHit  ctrl-c  to  quit\ni  =  %f\nt  =  %f\nFa  = 
%f \n  '  ,  i,  t.  Fa) 

pause 

end 

%This  checks  for  negative  tension: 
if  T (i, t)  <  0 
if  i  ==  1 

Told  =  T (i, t-1) ; 

else 

Told  =  T ( i-1 , t) ; 

end 

fprintf (' Tension  is  negative ! \nTry  making  less  drastic 
perturbations . \nNote :  the  code  cannot  analyze  slack 
conditions . \ni=%f\nt=%f\nT=%f\nT (previous ) =%f\nHit  ctrl-c  to 
quit\n' ,i,t,T(i,t) , Told) 
pause 

end 

%This  sets  values  for  the  characteristic  equations: 

Ga  =  W  (i,  t)  *sin  (Theta  (i,  t)  )  -  V  (i,  t)  *cos  (Theta  (i,  t)  )  + 
Fa*cos (Theta (i, t) ) ; 
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Gb  =  W  (i,  t)  *sin  (Theta  (i,  t)  )  -  V  (i,  t)  *cos  (Theta  (i,  t)  )  + 
Fb*cos (Theta (i, t) ) ; 

if  i  ==  L/ds+1  %Sets  values  at  body  -  based  on  body  mass  - 
this  is  currently  not  being  used 

Hb  =  (1/Mb)  *  (-FX  (i,  t)  -  WB*cos  (Phi  (i,  t)  )  )  ; 

Lb  =  (1/Mb)  *  (-FZ  (i,  t)  - 
WB*sin (Phi (i, t) ) *sin (Theta (i, t) ) ) ; 

else  %Sets  values  along  towline  -  based  on  towline  mass 
Hb  =  (1/mu) * (-FX (i, t)  -  Wt*cos (Phi (i, t) ) ) ; 

Lb  =  (1/mu)  *  (-FZ  (i,  t)  - 
Wt*sin (Phi (i, t) ) *sin (Theta (i, t) ) ) ; 
end 

Ha  =  Hb; 

Hatrack (i, t) =Ha; 

La  =  Lb; 

Ja  =  (V (i, t) -Fa) ; 

Jb  =  (V  (i,  t)  -Fb)  ; 

Ka  =  -U (i, t) *sin (Theta (i, t) ) ; 

Kb  =  Ka; 


%This  uses  interpolation  to  find  values  along  the  alpha 
characteristic : 

if  i  ~=  L/ds+1  %these  values  don't  exist  for  i=L/ds+l 

b  =  i*ds;  %b  position  is  the  next  position  along  the 


line 


a  =  (i-l)*ds;  %a  position  is  the  current  position  along 
the  line  (note:  line  STARTS  at  i=l,  thus  subtract  1) 

p  =  a  +  dt*Fa;  %comes  from  r-p  =  dt*Fa  and  the  fact 
that  r  =  a  (position  along  line) 

%  This  defines  values  for  interpolation: 

Ua  =  U(i,t);  %Ua  is  value  at  same  position  along  line 
w/  previous  timestep 

Ub  =  U(i+l,t);  %Ub  similarly  corresponds  to  the  next 
position  along  the  line 

Wa  =  W ( i , t ) ; 

Wb  =  W (i+1, t) ; 

Phia  =  Phi (i, t) ; 

Phib  =  Phi (i+1, t); 

Thetaa  =  Theta (i,t); 

Thetab  =  Theta (i+1, t) ; 

%  The  following  uses  straight  line  interpolation  to 

find  values 

Up  =  ( 1/ds) * ( (b-p) *Ua  +  (p-a) *Ub) ; 

Wp  =  ( 1/ds) *( (b-p) *Wa  +  (p-a) *Wb) ; 

Phip  =  (1/ds) *( (b-p) *Phia  +  (p-a) *Phib) ; 

Thetap  =  ( 1 /ds ) * ( (b-p) *Thetaa  +  (p-a) *Thetab) ; 

end 


%This  uses 
characteristic : 

if  i  ~ 


along  the  line 
the  line 


d 

a 


interpolation  to  find  values  along  the  beta 

1  %these  values  don't  exist  for  i=l 
(i-2)*ds;  %d  position  is  the  previous  position 

(i-l)*ds;  %a  position  is  the  current  position  along 


111 


q  =  a  +  dt*Fb;  %comes  from  r-q  =  dt*Fb  and  the  fact 
that  r  =  a  (position  along  line) 

%  This  defines  values  for  interpolation: 

Ua  =  U(i,t);  %Ua  is  value  at  same  position  along  line 
w/  previous  timestep 

Ud  =  U(i-l,t);  %Ub  similarly  corresponds  to  the 
previous  position  along  the  line 
Wa  =  W ( i , t ) ; 

Wd  =  W (i-1, t) ; 

Phia  =  Phi (i, t) ; 

Phid  =  Phi (i-1, t); 

Thetaa  =  Theta (i,t); 

Thetad  =  Theta (i-1, t) ; 

%  The  following  uses  straight  line  interpolation  to 

find  values 

Uq  =  (1/ds) * ( (q-d) *Ua  +  (a-q) *Ud) ; 

Wq  =  (1/ds) *( (q-d) *Wa  +  (a-q) *Wd) ; 

Phiq  =  (1/ds) *( (q-d) *Phia  +  (a-q) *Phid) ; 

Thetaq  =  (1/ds) *( (q-d) *Thetaa  +  (a-q) *Thetad) ; 

end 

%  This  finds  the  upper  limit  values  of  U,V,W,Phi,&  Theta: 
if  i  ==  1 

%  The  following  finds  phi  by  finding  a  zero  value  for 

the  fnct: 

phi  =  fzero(@(phi)  phi  -  (Up  +  Ga*Phip  -  Ha*dt  -  (C 
+  ustar) *sin (phi)  +  (D  +  vstar ) *cos (phi )) /Ga,  .001); 

%  The  following  finds  theta  by  finding  a  zero  value  for 

the  fnct: 

theta  =  f zero (@ (theta)  theta  -  (Wp  +  (C  + 
ustar) *sin (theta) *cos (phi)  -  (E  +  wstar) *cos (theta)  +  (D  + 
vstar) *sin (phi) *sin (theta)  +  Ja*Thetap  -  Ka* (phi  -  Phip)  -  La*dt)/Ja, 
.001)  ; 

if  theta  <  le-10  %due  to  round  off,  this  sets  to 

correct  zero  value 

theta  =  0; 

end 

Phi(i,t+1)  =  phi; 

Theta (i,t+l)  =  theta; 

U(i,t+1)  =  (C  +  ustar) *sin (phi)  -  (D  +  vstar ) *cos (phi ) ; 
V(i,t+1)  =  (C  +  ustar) *cos (theta) *cos (phi)  +  (D  + 
vstar) *cos (theta) *sin (phi)  +  (E  +  wstar) *sin (theta) ; 

W(i,t+1)  =  - (C  +  ustar) *sin (theta) *cos (phi)  -  (D  + 
vstar) *sin (phi) *sin (theta)  +  (E  +  wstar) *cos (theta) ; 
end 

%This  finds  the  values  down  the  entire  line  using  the 
characteristic  equations: 

if  i  ~=  1  &  i  ~=  L/ds+1 

Phi(i,t+1)  =  (Up  -  Uq  -  Gb*Phiq  +  Ga*Phip) / (Ga  -  Gb) ; 

U  (i, t+1)  =  (Up*Gb  -  Uq*Ga  +  Gb*Ga*(Phip  -  Phiq)  + 

Ha*dt* (Ga  -  Gb) ) / (Gb  -  Ga) ; 

W(i,t+1)  =  (Jb*Wp  -  Ja*Wq  +  Ja*Jb* (Thetap  -  Thetaq)  + 
Jb*Ka*(Phip  -  Phi (i,  t+1))  +  Ja*Kb* ( Phi ( i ,  t+1 )  -  Phiq)  +  La*dt* (Jb  - 
Ja)  )  /  (Jb  -  Ja)  ; 
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Theta (i,t+l)  =  (Thetap*Ja  -  Thetaq*Jb  +  Wp  -  Wq  - 
Phi (i, t+1 ) * (Ka  -  Kb)  +  Phip*Ka  -  Phiq*Kb) / (Ja  -  Jb)  ; 

V (i,  t+1)  =  V(i-l,t+l)  +  . 5* (U (i-1, t+1)  + 

U(i, t+1) ) *cos ( .5* (Theta (i-1, t+1)  +  Theta (i, t+1) ))* (Phi (i-1, t+1)  - 
Phi(i,t+1))  -  . 5* (W (i-1, t+1)  +  W (i, t+1) )* (Theta (i-1, t+1)  - 
Theta ( i ,  t+1 ) )  ; 

end 


%This  sets  up  the  aerodynamic  forces  on  the  line: 
if  i  ~=  L/ds+1  %Last  value  to  be  added  later 

F  =  Towf orce (U (i, t+1 ) , V (i, t+1 ) , W (i, t+1 ) , rho (i, t+1 ) ) ; 
FX ( i , t+1 )  =  - F  ( 1 ) ; 

FY ( i , t+1 )  =  -F (2 ) ; 

FZ (i, t+1)  =  -F  (3) ; 

end 

end 


%This  finds  the  values  at  lower  boundary  conditions: 


%Use  combination  of  equations  to  step  values  to  lower  boundary: 
Cl  =  Gb+V(L/ds, t+1) *cos (Theta (L/ds,  t+1)  )  - 
W(L/ds, t+1) *sin (Theta (L/ds,  t+1) )  +  (2*ds/dt)*cos (Theta (L/ds,  t+1) ) ; 
%constant  for  PhiL  denominator 
PhiL  = 

(  (2*ds/dt) *Phi (L/ds  +  1, t) *cos (Theta (L/ds, t+1) )  +  (V(L/ds,t+l) *cos (Theta ( L / 
ds, t+1) ) -W (L/ds, t+1) *sin (Theta (L/ds, t+1) ) )*(Phi(L/ds-l,t+l) ) -U (L/ds- 
1, t+1) -Hb*dt+Uq+Gb*Phiq) /Cl; 

ThetaL  =  ( (2*ds/dt) *Theta (L/ds+1, t) +V(L/ds, t+1) *Theta (L/ds- 
1, t+1) -U (L/ds, t+1) * (PhiL- Phi (L/ds-l,t+l) ) *sin (Theta (L/ds, t+1) ) -Wq- 
Jb*Thetaq+Lb*dt+W (L/ds-1, t+1) -Kb* (PhiL-Phiq) ) / (2*ds/dt+V (L/ds, t+1) -Jb) ; 
UL  =  Uq-Gb* (PhiL-Phiq) -Hb*dt; 

WL  =  Wq-Jb* (ThetaL-Thetaq) -Kb* (PhiL-Phiq) -Lb*dt; 

VL  =  V (L/ds-1, t+1) +W (L/ds, t+1) * (ThetaL- Theta (L/ds-1, t+1) ) - 
U(L/ds,t+l)*cos (Theta (L/ds, t+1) ) * (PhiL- Phi (L/ds-1, t+1) ) ; 

%Set  values: 

U (L/ds+1, t+1)  =  UL; 

V(L/ds+l, t+1)  =  VL; 

W (L/ds+1, t+1)  =  WL; 

Theta (L/ds+1 , t+1 )  =  ThetaL; 

Phi (L/ds+1 , t+1 )  =  PhiL; 


%Convert  to  Spatial: 

Spatial  vel  =  [UL  VL  WL] *Amatrix ( PhiL, ThetaL) ; 
uL  =  Spatial_vel ( 1 ) ; 
vL  =  Spatial  vel (2); 
wL  =  Spatial_vel(3); 
u (L/ds+1 , t+1 )  =  uL; 
v (L/ds+1, t+1)  =  vL; 
w (L/ds+1 , t+1 )  =  wL; 


%Find  drag  forces  on  the  body  and  tension  at  lower  end: 

F  =  Bodyforce (uL, vL, wL, rho (L/ds+1 , t+1 )) ;  %Set  forces  on 
body  in  spatial  axes 

if  t  ==  1  %Solving  for  t=2 : 

au  =  (u (L/ds+1, t+1) -u (L/ds+1, t) ) /dt; 


113 


av  =  (v (L/ds+1 , t+1 ) -v (L/ds+1 , t) ) /dt; 
aw  =  (w  (L/ds  +  1,  t+1) -w (L/ds  +  1, t) ) /dt; 

else 

au  =  (3*u (L/ds+1, t+1) -4*u (L/ds+1, t) +u (L/ds+1, t- 

1) ) / (2*dt)  ; 

av  =  (3*v (L/ds+1, t+1) -4*v (L/ds+1, t) +v (L/ds+1, t- 

1) ) / (2*dt) ; 

aw  =  (3*w (L/ds+1, t+1) -4*w (L/ds+1, t) +w (L/ds+1, t- 

1) ) / (2*dt) ; 

end 

Fx  =  -F(l)-Mb*au;  %negative  sign  on  force  is  due  to  drag 
force  orientation  being  same  direction  as  velocity 
Fy  =  -F (2) -WB-Mb*av; 

Fz  =  -F(3)-Mb*aw; 

FTL  =  Amatrix (PhiL, ThetaL) * [Fx; Fy; Fz] ; 

T (L/ds  +  1, t+1)  =  -FTL  (2) ; 

%Set  aerodynamic  forces  on  line  at  lower  end: 

F  =  Amatrix ( PhiL, ThetaL) * [F ( 1 ); F (2 ); F ( 3 )] ;  %Change  to 
towline  coordinates 

FX (L/ds+1, t+1)  =  - F ( 1 ) ; 

FY (L/ds+1, t+1)  =  -F (2 ) ; 

FZ (L/ds+1, t+1)  =  -F (3) ; 

%Set  tension  one  step  up  line: 

T(L/ds,t+l)  =  T (L/ds+1, t+1) +ds* ( (mu/dt) * (V(L/ds, t+1) -V(L/ds, t) - 
W(L/ds,  t+1) * (Theta (L/ds, t+1) -Theta (L/ds, t) ) +U (L/ds, t+1) 

*cos (Theta (L/ds, t+1) ) * (Phi (L/ds, t+1) -Phi (L/ds,  t) ) ) - 
FY (L/ds, t+1) +Wt*sin (Phi (L/ds, t+1) ) *cos (Theta (L/ds, t+1) ) ) ; 

%Set  tension  up  rest  of  line: 
for  i=L/ds-l : -1 : 1 

T (i, t+1)  =  T (i+2,  t+1) +2*ds* ( (mu/dt) * (V (i+1, t+1) -V (i+1, t) - 
W(i+l,t+l)* (Theta (i+1, t+1) -Theta (i+1, t) ) +U (i+1, t+1) 

*cos (Theta (i  +  1,  t+1) ) * (Phi (i  +  1,  t+1) -Phi (i+1,  t) )  )  - 
FY (i+1, t+1) +Wt*sin (Phi (i+1, t+1) ) *cos (Theta (i+1, t+1) ) ) ; 
end 

for  i=l:l:L/ds+l 

%Find  spatial  velocity  values: 

Spatial_vel  =  [U(i,t+1)  V(i,t+1) 

W ( i , t+1 ) ] * Amatrix ( Phi ( i , t+1 )  ,  Theta ( i ,  t+1 )  )  ; 

u(i,t+l)  =  Spatial_vel (1) ; 
v(i,t+l)  =  Spatial_vel(2); 
w(i,t+l)  =  Spatial_vel ( 3 ) ; 

%Find  spatial  position  values: 

if  i  ==  1  %Sets  the  value  at  the  a/c  -  subtracing  the 
initial  velocities  keeps  system  about  these  velocities 
x(i,t+l)  =  x(i,t)  +  (u (i, t+1) -C) *dt; 
y ( i , t+1 )  =  y (i, t)  +  (v (i, t+1) -D) *dt; 
z(i,t+l)  =  z(i,t)  +  (w (i, t+1) -E) *dt; 

else 

x(i,t+l)  =  x(i-l,t+l)  - 
ds*cos (Theta (i, t+1) ) *cos (Phi (i, t+1) ) ; 
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y (i, t+1)  =  y (i-1, t+1)  - 
ds*cos (Theta (i, t+1) ) *sin (Phi (i, t+1) ) ; 

z(i,t+l)  =  z(i-l,t+l)  -  ds*sin  (Theta (i, t+1) ) ; 

end 

%The  following  is  used  to  find  heat  transfer  and  density  in 
ambient  temp. 

%Due  to  computing  power,  it  can  be  commented  out  when  not 
analyzing  heat  transfer. 

Tempset; 


end 

end 
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Appendix  C:  MATLAB®  Code  -  Tempset  Code 


%This  takes  in  values  of  x,  y,  and  z  in  meters,  U,  V,  and  W  as  m/s, 
%Told  in  K,  and  t  as  timestep  and  outputs  the  temperature  and  air 
%density  at  the  current  position  for  the  next  timestep.  This  method 
%can  be  applied  to  the  entire  towline.  It  is  not  applied  to  the  towed 
%body,  but  values  are  calculated  at  the  point  of  attachment  between  the 
%line  and  body,  and  these  values  are  assumed  to  be  the  density  values 
%for  the  body.  The  body  should  be  sufficiently  far  from  a  heat  source, 
%thus  this  method  should  give  accurate  data. 

o, 

o 

%Time  spacing  in  Temp^set  is  same  as  dt,  thus  ' t '  is  timestep  value. 

o, 

o 

%x,y,z  values  are  relative  to  the  initial  aircraft  position  (this  can 
be 

%altered  by  setting  x,y,z  values  in  other  code  to  be  dependent  on  a 
%different  position) . 

%Set  up  coordinates  for  interpolation: 

%The  +1  value  is  due  to  the  first  grid  position  being  at  zero 
xpos  =  x (i, t+1) *GS+l-xoffset*GS;  %sets  true  x  position  within  grid 
points 

ypos  =  y (i, t+1) *GS+l-yoffset*GS;  %sets  true  y  position  within  grid 
points 

zpos  =  z (i, t+1) *GS+l-zoffset*GS;  %sets  true  z  position  within  grid 
points 

xl  =  floor (xpos);  %finds  the  lower  position  for  x  for  Temp  set 
x2  =  ceil (xpos);  %finds  the  upper  position  for  x  for  Temp  set 
yl  =  floor (ypos);  %finds  the  lower  position  for  y  for  Temp  set 
y2  =  ceil (ypos);  %finds  the  upper  position  for  y  for  Temp_set 
zl  =  floor(zpos);  %finds  the  lower  position  for  z  for  Temp  set 
z2  =  ceil (zpos);  %finds  the  upper  position  for  z  for  Temp  set 

%The  'if'  statements  prevent  an  output  of  a  zero  value  later  due  to 
xl=x2,  yl=y2,  or  zl=z2. 
if  xl==x2 

x2=xl+l ; 

end 

if  yl==y2; 

y2=yl+l ; 

end 

if  zl==z2; 

z2=zl+l; 

end 

%Set  up  four  values  for  interpolation  at  eight  points: 

T1  =  Temp_set (xl, yl, zl, t+1) ; 

T2  =  Temp  set (x2, yl, zl, t+1) ; 

T3  =  Temp_set (xl , y2 , zl , t+1 ) ; 

T4  =  Temp  set (x2, y2, zl, t+1) ; 

T5  =  Temp_set (xl , yl , z2 , t+1 ) ; 

T6  =  Temp_set (x2 , yl , z2 , t+1 ) ; 
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T7  =  Temp  set (xl , y2 , z2 , t+1 ) ; 
T8  =  Temp_set (x2 , y2 , z2 , t+1 ) ; 


%Set  values  at  current  position  through  interpolation: 

%Since  everything  is  in  terms  of  grid  position  (i.e.,  not  real 
position),  no  need  to  divide  by  grid  size. 

Txl  =  (x2-xpos)*Tl  +  (xpos-xl ) *T2 ;  %value  along  x-axis  at  yl,zl 
Tx2  =  (x2-xpos)*T3  +  (xpos-xl ) *T4 ;  %value  along  x-axis  at  y2,zl 
Tx3  =  (x2-xpos)*T5  +  (xpos-xl ) *T6;  %value  alone  x-axis  at  yl,z2 
Tx4  =  (x2-xpos)*T7  +  (xpos-xl ) *T8 ;  %value  alone  x-axis  at  y2,z2 
Txyl  =  (y2-ypos) *Txl  +  (ypos-yl ) *Tx2 ;  %value  in  x-y  plane  at  zl 
Txy2  =  (y2-ypos) *Tx3  +  (ypos-yl ) *Tx4 ;  %value  in  x-y  plane  at  z2 
Txyz  =  (z2-zpos) *Txyl  +  ( zpos-zl ) *Txy2 ;  %final  air  temp 

interpolated  at  the  position 


%This  should  be  commented  out  for  now  since  the  Nusselt  number  in 
%currently  only  for  perpendicular  velocity  (i.e.,  heat  transfer  only 
%happens  in  perpendicular  direction) : 


%Set  the  diameter  to  calculate  values  (we  only  use  perpendicular 
velocity  due  to  Nusselt  number  calculations  -  this  is  left  for  future 
work) : 

%  Vperp  =  sqrt (U ( i , t+1 ) A2  +  W(i,t+1)A2);  %perpendicular  velocity 

%  Vpar  =  abs (V ( i , t+1 ) ) ;  %parallel  velocity  -  make  positiv 

%  if  Vpar  ~=  0;  %divide  by  zero  error  if  it  does 

%  Angle  =  atan (Vperp/Vpar ) ;  %angle  at  which  airflow  acts 

%  if  Vperp/Vpar  >  params . dL/set . ds  %if  the  air  completely 

crosses  the  perpendicular  component  of  line  over  the  interval  ds 
%  DL  =  params . dL/sin (Angle) ;  %diameter  over  which  airflow 

acts 

%  else  %if  the  air  completely  crosses  the  parallel  component  of 

line  over  the  interval  ds 

%  DL  =  set . ds/cos (Angle) ;  %diameter  over  which  airflow  acts 

%  end 

%  else 

%  DL  =  set.ds;  %diameter  over  which  airflow  acts 

%  end 


DL  =  params.dL; 
diameter 


^diameter  over  which  airflow  acts  -  perpendicular 


%Calculate  values  to  find  temp  and  rho : 

%01d  values: 

Tf= (Temp ( i , t+1 ) +set . CT) /2 ;  %film  temp  based  on  old  line  temp 
nu  =  8e-10* (Tf ) A1 . 7235;  %kinematic  viscosity  based  on  old  temp 
%Calculate  values  based  on  old  temp: 

Re  =  abs (U (i, t+1) ) *DL/nu;  %Reynolds  number  -  this  would  change 
for  a  different  Nu  number  calculation 

Nu  =  0.42*PrA0.2  +  . 057*PrA (1/3) *sqrt (Re) ;  %Nusselt  number 
k  =  0 . 0002235* (Tf) A0 . 8302 ;  %thermal  conductivity  of  air 
h  =  Nu*k/params .dL;  %convection  coefficient 
%Calculate  new  values: 

Temp(i,t+2)  =  Temp (i, t+1)  +  Area*h*dt* (Tf- 
Temp ( i , t+1 ) ) / (Vol*params . rhoL*Cp) ;  %new  temp  on  line 
Tf  =  (Temp ( i ,  t+2 ) +set . CT) /2 ;  %new  film  temp 

rho(i,t+2)  =  357 . 88* (Tf ) A (-1 . 0041) ;  %new  air  density  in  film 
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Appendix  D:  MATLAB®  Code  -  Other  Functions 


The  first  function  is  used  to  set  the  initial  values  for  the  entire  simulation,  and  is 
called  up  in  multiple  other  functions.  It  is  currently  set  to  model  a  1  second  vertical 
acceleration  to  30  m/s,  with  an  immediate  acceleration  back  to  zero  over  1  second. 

function  set=set_values ( ) 

%Set  initial  conditions: 

set.u  =  100;  %intial  vel  in  x-dir 
set.v  =  0;  %initial  vel  in  y-dir 
set.w  =  0;  %initial  vel  in  z-dir 

set. Phi  =  0.01;  %initial  angle  of  towline  (constant  down  line)  - 
nonzero  prevents  negative  tension 

set. Theta  =  0;  %initial  angle  of  towline  (constant  down  line) 

%Set  mesh  spacing: 

set.ds  =  5;  %change  in  length  (meters) 
set.dt  =  .01;  %change  in  time  (sec) 
set. time  =  10;  %sec  (total  time) 

%First  perturbation: 

set.upert  =  0.0;  %perturbations  in  velocity  of  a/c 
set.vpert  =  30.0; 

set.wpert  =  0.0;  %limits  in  these  values 

set.tpert  =  1;  ^constant  acceleration  until  reaches  final 
perturbation  value  at  tpert 

%Second  perturbnation : 

set.tstart2  =  1;  %time  to  start  second  perturbation 

set.upert2  =  0.0; 

set.vpert2  =  -30.0; 

set.wpert2  =  0.0; 

set.tpert2  =  1; 

%Set  constant  temp  value  [K]  for  air  for  Tempset.m  (to  be  changed  for 
temp  field  later) 

set.CT  =  288.2;  %Sea  level  St.  Atm.  and  Press,  value  is  288.2  K 

%Set  value  for  initial  line  temp  [K] : 
set.LTemp  =  288.2; 


The  second  function  was  modified  from  Richardson  and  sets  the  initial  parameters 
along  the  line.  These  values  were  set  somewhat  arbitrarily  (although  some  values,  such 
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as  the  density  of  steel,  were  set  specifically)  to  analyze  different  towline  behaviors.  This 
function  is  also  called  into  Richardson’s  work  for  steady  state  analysis: 

function  params=line  params ( ) 

%Original  program  written  by  Ralph  Anthenien  and  Tyler  Richardson, 
AFIT/ENY 

%Modified  by  Ralph  Anthenien  and  Christopher  Hill,  AFIT/ENY 
%function  line_params 

%return  parameters  about  towline  and  decoy 
set=set_values; 

params . g=9 . 8 ;  %gravitational  acceleration  [m/sA2] 
params .dB=0 . 496;  %decoy  diam 

params . Ad=pi*params .dBA2/4;  %frontal  area  [mA2] 
params . CdD=l ;  %Decoy  drag  coeff 
params. mD=4;  %decoy  mass  [kg] 

params . WB=params . mD*params . g;  %decoy  weight  [N] 
params. LD=1;  %Decoy  length  [m] 

params . CdL=l . 1 ;  %Perpendicular  line  drag  coeff 

params . CdLSF=0 . 04 ;  %Skin  friction  line  drag  coeff  (for  Y-axis) 
params . rhoL=7600;  %line  density  [kg/mA3]  -  approx,  value  for  most 
steels 

params . dL=0 . 00127 ;  %Line  Diam  [m] 

params . mu=params . rhoL*pi*params . dLA2 /4 ;  %Mass  of  towline  per  meter 
length  [kg/m] 

params .Wt=params .mu*params . g;  %Weight  of  towline  per  meter  [N] 
params . Vx=set . u;  %velocity  (lookup  to  be  used  later)  [m/s] 
params .Vy=set .v;  %upward  vel  [m/s] 
params .Vz=set .w;  %transverse  vel  [m/s] 

params . V=sqrt (params . VxA2+params . VyA2+params . Vz A2 ) ;  %relative  vel 
params . vchar=params . V;  %characteristic  velocity  to  calc  char  drag 
tension  -  assumed  to  be  same  as  total  velocity 

params . rhoa=357 .88* ( ( set . LTemp+set .CT)/2)A (-1.0041);  %air  density 
[kg/mA3] 

params . LL=30 ;  %line  length  [m] 

params. tc=l;  %time  constant  (for  transient  only) 

params . T0=params . LL*params . g*pi*params . rhoL*params . dLA2 / 4 ; 


The  third  function  is  the  transformation  matrix  to  convert  from  towline  to  space 
coordinates  and  vice  versa: 

function  A  =  Amatrix (phi ,  theta) 

%  This  transforms  from  space  to  towline  coordinate  systems  as  follows: 
%  Towline (3x1)  =  [A] *Space (3x1) 

%  Space (3x1)  =  (Towline (1x3) * [A] ) ' 

%  note:  phi  and  theta  are  angles  btwn  current  and  old  line  positions 
All  =  sin (phi ) ; 

A12  =  -cos (phi ) ; 

A13  =  0; 

A21  =  cos (phi) *cos (theta); 
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A22  =  cos  (theta) *sin (phi) ; 

A23  =  sin (theta); 

A31  =  -cos (phi) *sin (theta) ; 

A32  =  -sin (theta) *sin (phi) ; 

A33  =  cos (theta) ; 

A  =  [All  A12  A13;  A21  A22  A23;  A31  A32  A33]; 


The  fourth  function  calculates  the  forces  on  the  body  in  spatial  coordinates: 

function  F  =  Bodyforce (u, v, w, rho) 

%  This  calculates  the  total  aerodynamic  forces  on  the  body 
params=line  params; 

Tot  vel  =  sqrt (uA2  +  vA 2  +  wA2); 

Dx  =  .  5*rho*u*Tot  vel*params . CdD*params . Ad;  %Drag  in  x-direction 

Dy  =  . 5*rho*v*Tot  vel*params . CdD*params . Ad;  %Drag  in  y-direction 

Dz  =  . 5*rho*w*Tot  vel*params . CdD*params . Ad;  %Drag  in  z-direction 

F= [Dx  Dy  Dz]  ; 


The  fifth  function  calculates  the  forces  on  the  towline  in  towline  coordinates: 

function  F  =  Towforce (U, V, W, rho) 

%  This  calculates  the  aerodynamic  forces  on  the  towline  per  meter 
params=line  params; 

V_perp  =  sqrt (UA2  +  WA2);  %total  velocity  in  the  perpendicular 
direction 

DX  =  . 5*rho*params . dL*U*V  perp* (params . CdL+params . CdLSF) ;  %Drag  in  X 
direction 

DY  =  0;  %Drag  in  Y-direction 

DZ  =  . 5*rho*params . dL*W*V  perp* (params . CdL+params . CdLSF) ;  %Drag  in  Z 

direction 

F= [DX  DY  DZ] ; 
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